NDFS Synthesis Manuscript: Chlorophyll analysis
Model selection
Dave Bosworth
August 31, 2023
Purpose
Explore and analyze the continuous chlorophyll data to be included in the NDFS synthesis manuscript. We will attempt to fit multiple models to the data set and use a model selection process to determine the best one.
Global code and functions
# Load packages
library(tidyverse)
library(scales)
library(knitr)
library(mgcv)
library(car)
library(gratia)
library(here)
library(conflicted)
# Source functions
source(here("manuscript_synthesis/src/global_functions.R"))
# Declare package conflict preferences
conflicts_prefer(dplyr::filter())Display current versions of R and packages used for this analysis:
devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.2.3 (2023-03-15 ucrt)
## os Windows 10 x64 (build 19044)
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate English_United States.utf8
## ctype English_United States.utf8
## tz America/Los_Angeles
## date 2023-08-31
## pandoc 3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 4.2.0)
## bslib 0.4.2 2022-12-16 [1] CRAN (R 4.2.2)
## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.2.3)
## callr 3.7.3 2022-11-02 [1] CRAN (R 4.2.2)
## car * 3.1-2 2023-03-30 [1] CRAN (R 4.2.3)
## carData * 3.0-5 2022-01-06 [1] CRAN (R 4.2.1)
## cli 3.6.1 2023-03-23 [1] CRAN (R 4.2.3)
## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.2.2)
## conflicted * 1.2.0 2023-02-01 [1] CRAN (R 4.2.2)
## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.2.1)
## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.2.1)
## digest 0.6.31 2022-12-11 [1] CRAN (R 4.2.2)
## dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.2.3)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.1)
## evaluate 0.21 2023-05-05 [1] CRAN (R 4.2.3)
## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.2.2)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.2.2)
## forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.2.2)
## fs * 1.6.2 2023-04-25 [1] CRAN (R 4.2.3)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.1)
## ggplot2 * 3.4.2 2023-04-03 [1] CRAN (R 4.2.3)
## glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.1)
## gratia * 0.8.1 2023-02-02 [1] CRAN (R 4.2.2)
## gtable 0.3.3 2023-03-21 [1] CRAN (R 4.2.3)
## here * 1.0.1 2020-12-13 [1] CRAN (R 4.2.1)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.2.3)
## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.2.3)
## htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.2.3)
## httpuv 1.6.9 2023-02-14 [1] CRAN (R 4.2.2)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.2.1)
## jsonlite 1.8.4 2022-12-06 [1] CRAN (R 4.2.2)
## knitr * 1.42 2023-01-25 [1] CRAN (R 4.2.2)
## later 1.3.0 2021-08-18 [1] CRAN (R 4.2.1)
## lattice 0.21-8 2023-04-05 [1] CRAN (R 4.2.3)
## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.2.1)
## lubridate * 1.9.2 2023-02-10 [1] CRAN (R 4.2.2)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.1)
## Matrix 1.5-4 2023-04-04 [1] CRAN (R 4.2.3)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.2.1)
## mgcv * 1.8-42 2023-03-02 [1] CRAN (R 4.2.2)
## mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.2.1)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.1)
## mvnfast 0.2.8 2023-02-23 [1] CRAN (R 4.2.2)
## nlme * 3.1-162 2023-01-31 [2] CRAN (R 4.2.3)
## patchwork 1.1.2 2022-08-19 [1] CRAN (R 4.2.1)
## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.3)
## pkgbuild 1.4.0 2022-11-27 [1] CRAN (R 4.2.2)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.1)
## pkgload 1.3.2 2022-11-16 [1] CRAN (R 4.2.2)
## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.2.1)
## processx 3.8.1 2023-04-18 [1] CRAN (R 4.2.3)
## profvis 0.3.7 2020-11-02 [1] CRAN (R 4.2.1)
## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.2.1)
## ps 1.7.5 2023-04-18 [1] CRAN (R 4.2.3)
## purrr * 1.0.1 2023-01-10 [1] CRAN (R 4.2.2)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.1)
## Rcpp 1.0.10 2023-01-22 [1] CRAN (R 4.2.2)
## readr * 2.1.4 2023-02-10 [1] CRAN (R 4.2.2)
## remotes 2.4.2 2021-11-30 [1] CRAN (R 4.2.1)
## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.2.3)
## rmarkdown 2.21 2023-03-26 [1] CRAN (R 4.2.3)
## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.2.1)
## rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.1)
## sass 0.4.6 2023-05-03 [1] CRAN (R 4.2.3)
## scales * 1.2.1 2022-08-20 [1] CRAN (R 4.2.1)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.1)
## shiny 1.7.4 2022-12-15 [1] CRAN (R 4.2.2)
## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.2.2)
## stringr * 1.5.0 2022-12-02 [1] CRAN (R 4.2.2)
## tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.2.3)
## tidyr * 1.3.0 2023-01-24 [1] CRAN (R 4.2.2)
## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.1)
## tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.2.2)
## timechange 0.2.0 2023-01-11 [1] CRAN (R 4.2.2)
## tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.2.3)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.2.1)
## usethis 2.1.6 2022-05-25 [1] CRAN (R 4.2.1)
## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.2.2)
## vctrs 0.6.2 2023-04-19 [1] CRAN (R 4.2.3)
## withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.1)
## xfun 0.39 2023-04-20 [1] CRAN (R 4.2.3)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.2.1)
## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.2.2)
##
## [1] C:/R/win-library/4.2
## [2] C:/Program Files/R/R-4.2.3/library
##
## ──────────────────────────────────────────────────────────────────────────────Import Data
# Define file path for processed data
fp_data <- here("manuscript_synthesis/data/processed")
# Import daily average chlorophyll and water temperature data
df_chla_wt <- readRDS(file.path(fp_data, "chla_wt_daily_avg_2013-2019.rds"))
# Import daily average flow data
df_flow <- readRDS(file.path(fp_data, "flow_daily_avg_2013-2019.rds"))Prepare Data
# Create a vector for the factor order of StationCode
sta_order <- c(
"RCS",
"RD22",
"I80",
"LIS",
"STTD",
"LIB",
"RVB"
)
# We will use LIS flow data as a proxy for STTD and RD22 flow data as a proxy
# for I-80
df_flow_c <- df_flow %>%
filter(StationCode %in% c("RD22", "LIS")) %>%
mutate(StationCode = case_match(StationCode, "RD22" ~ "I80", "LIS" ~ "STTD")) %>%
bind_rows(df_flow)
# Prepare chlorophyll and flow data for exploration and analysis
df_chla_wt_c1 <- df_chla_wt %>%
# Join flow data to chlorophyll data
left_join(df_flow_c, by = join_by(StationCode, Date)) %>%
mutate(
# Scale and log transform chlorophyll values
Chla_log = log(Chla * 1000),
# Add Year and DOY variables
Year = year(Date),
DOY = yday(Date)
) %>%
# Add Flow Action Periods
ndfa_action_periods() %>%
mutate(
# Apply factor orders to FlowActionPeriod and StationCode
FlowActionPeriod = factor(FlowActionPeriod, levels = c("Before", "During", "After")),
StationCode = factor(StationCode, levels = sta_order),
# Add a column for Year as a factor for the model
Year_fct = factor(Year)
) %>%
arrange(StationCode, Date)Explore sample counts by Station
df_chla_wt_c1 %>%
summarize(
min_date = min(Date),
max_date = max(Date),
num_samples = n(),
.by = c(StationCode, Year, FlowActionPeriod)
) %>%
arrange(StationCode, Year, FlowActionPeriod) %>%
kable()| StationCode | Year | FlowActionPeriod | min_date | max_date | num_samples |
|---|---|---|---|---|---|
| RCS | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| RCS | 2014 | During | 2014-09-09 | 2014-09-23 | 11 |
| RCS | 2014 | After | 2014-09-24 | 2014-09-24 | 1 |
| RCS | 2015 | Before | 2015-08-17 | 2015-08-20 | 4 |
| RCS | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RCS | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| RCS | 2016 | Before | 2016-06-23 | 2016-07-13 | 20 |
| RCS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RCS | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RCS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RCS | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RCS | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RCS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RCS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RCS | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RCS | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RCS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RCS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| RD22 | 2013 | Before | 2013-08-15 | 2013-08-21 | 7 |
| RD22 | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| RD22 | 2013 | After | 2013-10-03 | 2013-11-04 | 33 |
| RD22 | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| RD22 | 2014 | During | 2014-09-09 | 2014-09-23 | 15 |
| RD22 | 2014 | After | 2014-09-24 | 2014-11-08 | 46 |
| RD22 | 2015 | Before | 2015-07-24 | 2015-08-20 | 28 |
| RD22 | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RD22 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| RD22 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| RD22 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RD22 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RD22 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RD22 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RD22 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RD22 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RD22 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RD22 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RD22 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RD22 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RD22 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| I80 | 2013 | Before | 2013-08-15 | 2013-08-21 | 7 |
| I80 | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| I80 | 2013 | After | 2013-10-03 | 2013-11-01 | 30 |
| I80 | 2014 | Before | 2014-07-25 | 2014-09-08 | 43 |
| I80 | 2014 | During | 2014-09-09 | 2014-09-23 | 15 |
| I80 | 2014 | After | 2014-09-24 | 2014-11-08 | 46 |
| I80 | 2015 | Before | 2015-07-24 | 2015-08-07 | 15 |
| I80 | 2015 | During | 2015-08-27 | 2015-10-01 | 36 |
| I80 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| I80 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| I80 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| I80 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| I80 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| I80 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| I80 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| I80 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| I80 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| I80 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| I80 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| I80 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| I80 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIS | 2013 | Before | 2013-07-07 | 2013-08-21 | 39 |
| LIS | 2013 | During | 2013-08-22 | 2013-10-02 | 38 |
| LIS | 2013 | After | 2013-10-03 | 2013-11-17 | 34 |
| LIS | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| LIS | 2014 | During | 2014-09-09 | 2014-09-23 | 14 |
| LIS | 2014 | After | 2014-09-24 | 2014-11-08 | 46 |
| LIS | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIS | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| LIS | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIS | 2016 | Before | 2016-05-29 | 2016-07-13 | 46 |
| LIS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIS | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| LIS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIS | 2017 | During | 2017-08-29 | 2017-09-09 | 12 |
| LIS | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| LIS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| LIS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| LIS | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| LIS | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| LIS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| LIS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| STTD | 2013 | Before | 2013-08-15 | 2013-08-21 | 7 |
| STTD | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| STTD | 2013 | After | 2013-10-03 | 2013-11-04 | 33 |
| STTD | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| STTD | 2014 | During | 2014-09-09 | 2014-09-23 | 15 |
| STTD | 2014 | After | 2014-09-24 | 2014-11-08 | 29 |
| STTD | 2015 | Before | 2015-07-27 | 2015-08-20 | 25 |
| STTD | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| STTD | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| STTD | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| STTD | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| STTD | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| STTD | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| STTD | 2017 | During | 2017-08-29 | 2017-09-01 | 4 |
| STTD | 2017 | After | 2017-09-20 | 2017-09-25 | 6 |
| STTD | 2018 | Before | 2018-07-20 | 2018-08-27 | 39 |
| STTD | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| STTD | 2018 | After | 2018-09-27 | 2018-10-15 | 19 |
| STTD | 2019 | Before | 2019-07-26 | 2019-08-25 | 31 |
| STTD | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| STTD | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIB | 2013 | Before | 2013-07-16 | 2013-08-21 | 37 |
| LIB | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| LIB | 2013 | After | 2013-10-03 | 2013-11-16 | 45 |
| LIB | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| LIB | 2014 | During | 2014-09-09 | 2014-09-23 | 15 |
| LIB | 2014 | After | 2014-09-24 | 2014-11-05 | 43 |
| LIB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIB | 2015 | During | 2015-08-21 | 2015-10-01 | 38 |
| LIB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIB | 2016 | Before | 2016-05-29 | 2016-07-13 | 46 |
| LIB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIB | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| LIB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| LIB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| LIB | 2018 | Before | 2018-08-14 | 2018-08-27 | 14 |
| LIB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| LIB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| LIB | 2019 | Before | 2019-07-11 | 2019-07-30 | 20 |
| LIB | 2019 | During | 2019-09-05 | 2019-09-21 | 17 |
| LIB | 2019 | After | 2019-09-22 | 2019-11-06 | 39 |
| RVB | 2013 | Before | 2013-07-07 | 2013-08-21 | 46 |
| RVB | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| RVB | 2013 | After | 2013-10-03 | 2013-11-17 | 46 |
| RVB | 2014 | Before | 2014-07-25 | 2014-09-07 | 45 |
| RVB | 2014 | During | 2014-09-10 | 2014-09-23 | 14 |
| RVB | 2014 | After | 2014-09-24 | 2014-11-08 | 46 |
| RVB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| RVB | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RVB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| RVB | 2016 | Before | 2016-05-29 | 2016-07-13 | 39 |
| RVB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RVB | 2016 | After | 2016-08-02 | 2016-09-16 | 37 |
| RVB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RVB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RVB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RVB | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RVB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RVB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RVB | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RVB | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RVB | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
Looking at the sample counts and date ranges, we’ll only include years 2015-2019 for the analysis. However, there are still a few Station-Year-FlowPeriod combinations with low sample counts within 2015-2019.
df_chla_wt_c2 <- df_chla_wt_c1 %>%
filter(Year %in% 2015:2019) %>%
mutate(Year_fct = fct_drop(Year_fct))
df_chla_wt_c2 %>%
summarize(
min_date = min(Date),
max_date = max(Date),
num_samples = n(),
.by = c(StationCode, Year, FlowActionPeriod)
) %>%
arrange(StationCode, Year, FlowActionPeriod) %>%
kable()| StationCode | Year | FlowActionPeriod | min_date | max_date | num_samples |
|---|---|---|---|---|---|
| RCS | 2015 | Before | 2015-08-17 | 2015-08-20 | 4 |
| RCS | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RCS | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| RCS | 2016 | Before | 2016-06-23 | 2016-07-13 | 20 |
| RCS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RCS | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RCS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RCS | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RCS | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RCS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RCS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RCS | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RCS | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RCS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RCS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| RD22 | 2015 | Before | 2015-07-24 | 2015-08-20 | 28 |
| RD22 | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RD22 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| RD22 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| RD22 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RD22 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RD22 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RD22 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RD22 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RD22 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RD22 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RD22 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RD22 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RD22 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RD22 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| I80 | 2015 | Before | 2015-07-24 | 2015-08-07 | 15 |
| I80 | 2015 | During | 2015-08-27 | 2015-10-01 | 36 |
| I80 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| I80 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| I80 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| I80 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| I80 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| I80 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| I80 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| I80 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| I80 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| I80 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| I80 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| I80 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| I80 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIS | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIS | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| LIS | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIS | 2016 | Before | 2016-05-29 | 2016-07-13 | 46 |
| LIS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIS | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| LIS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIS | 2017 | During | 2017-08-29 | 2017-09-09 | 12 |
| LIS | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| LIS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| LIS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| LIS | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| LIS | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| LIS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| LIS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| STTD | 2015 | Before | 2015-07-27 | 2015-08-20 | 25 |
| STTD | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| STTD | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| STTD | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| STTD | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| STTD | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| STTD | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| STTD | 2017 | During | 2017-08-29 | 2017-09-01 | 4 |
| STTD | 2017 | After | 2017-09-20 | 2017-09-25 | 6 |
| STTD | 2018 | Before | 2018-07-20 | 2018-08-27 | 39 |
| STTD | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| STTD | 2018 | After | 2018-09-27 | 2018-10-15 | 19 |
| STTD | 2019 | Before | 2019-07-26 | 2019-08-25 | 31 |
| STTD | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| STTD | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIB | 2015 | During | 2015-08-21 | 2015-10-01 | 38 |
| LIB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIB | 2016 | Before | 2016-05-29 | 2016-07-13 | 46 |
| LIB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIB | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| LIB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| LIB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| LIB | 2018 | Before | 2018-08-14 | 2018-08-27 | 14 |
| LIB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| LIB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| LIB | 2019 | Before | 2019-07-11 | 2019-07-30 | 20 |
| LIB | 2019 | During | 2019-09-05 | 2019-09-21 | 17 |
| LIB | 2019 | After | 2019-09-22 | 2019-11-06 | 39 |
| RVB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| RVB | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RVB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| RVB | 2016 | Before | 2016-05-29 | 2016-07-13 | 39 |
| RVB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RVB | 2016 | After | 2016-08-02 | 2016-09-16 | 37 |
| RVB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RVB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RVB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RVB | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RVB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RVB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RVB | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RVB | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RVB | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
Before we run our models, we’ll need to remove the NA
values in the Flow variable. When we remove all records where there are
missing daily average flow values, there are even more
Station-Year-FlowPeriod combinations with low sample counts within
2015-2019.
df_chla_wt_q <- df_chla_wt_c2 %>% drop_na(Flow)
df_chla_wt_q %>%
drop_na(Flow) %>%
summarize(
min_date = min(Date),
max_date = max(Date),
num_samples = n(),
.by = c(StationCode, Year, FlowActionPeriod)
) %>%
arrange(StationCode, Year, FlowActionPeriod) %>%
kable()| StationCode | Year | FlowActionPeriod | min_date | max_date | num_samples |
|---|---|---|---|---|---|
| RCS | 2015 | Before | 2015-08-17 | 2015-08-20 | 4 |
| RCS | 2015 | During | 2015-08-21 | 2015-09-28 | 39 |
| RCS | 2016 | Before | 2016-06-23 | 2016-07-13 | 8 |
| RCS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RCS | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RCS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RCS | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RCS | 2017 | After | 2017-09-19 | 2017-11-03 | 36 |
| RCS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RCS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RCS | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RCS | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RCS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RCS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| RD22 | 2015 | Before | 2015-07-24 | 2015-08-20 | 28 |
| RD22 | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RD22 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| RD22 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| RD22 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RD22 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RD22 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RD22 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RD22 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RD22 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RD22 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RD22 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RD22 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RD22 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RD22 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| I80 | 2015 | Before | 2015-07-24 | 2015-08-07 | 15 |
| I80 | 2015 | During | 2015-08-27 | 2015-10-01 | 36 |
| I80 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| I80 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| I80 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| I80 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| I80 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| I80 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| I80 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| I80 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| I80 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| I80 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| I80 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| I80 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| I80 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIS | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIS | 2015 | During | 2015-08-21 | 2015-10-01 | 37 |
| LIS | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIS | 2016 | Before | 2016-05-29 | 2016-07-13 | 46 |
| LIS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIS | 2016 | After | 2016-08-02 | 2016-09-16 | 43 |
| LIS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIS | 2017 | During | 2017-08-29 | 2017-09-09 | 12 |
| LIS | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| LIS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| LIS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| LIS | 2018 | After | 2018-09-27 | 2018-11-11 | 42 |
| LIS | 2019 | Before | 2019-07-19 | 2019-08-25 | 34 |
| LIS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| LIS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| STTD | 2015 | Before | 2015-07-27 | 2015-08-20 | 25 |
| STTD | 2015 | During | 2015-08-21 | 2015-10-01 | 37 |
| STTD | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| STTD | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| STTD | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| STTD | 2016 | After | 2016-08-02 | 2016-09-16 | 43 |
| STTD | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| STTD | 2017 | During | 2017-08-29 | 2017-09-01 | 4 |
| STTD | 2017 | After | 2017-09-20 | 2017-09-25 | 6 |
| STTD | 2018 | Before | 2018-07-20 | 2018-08-27 | 39 |
| STTD | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| STTD | 2018 | After | 2018-09-27 | 2018-10-11 | 15 |
| STTD | 2019 | Before | 2019-07-26 | 2019-08-25 | 27 |
| STTD | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| STTD | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIB | 2015 | During | 2015-08-21 | 2015-10-01 | 33 |
| LIB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIB | 2016 | Before | 2016-05-29 | 2016-07-13 | 36 |
| LIB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIB | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| LIB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| LIB | 2017 | After | 2017-09-19 | 2017-11-02 | 29 |
| LIB | 2018 | Before | 2018-08-14 | 2018-08-27 | 9 |
| LIB | 2018 | During | 2018-08-28 | 2018-09-08 | 12 |
| LIB | 2018 | After | 2018-09-27 | 2018-10-31 | 15 |
| LIB | 2019 | Before | 2019-07-11 | 2019-07-30 | 18 |
| LIB | 2019 | During | 2019-09-07 | 2019-09-21 | 10 |
| LIB | 2019 | After | 2019-09-22 | 2019-10-25 | 3 |
| RVB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| RVB | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RVB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| RVB | 2016 | Before | 2016-05-29 | 2016-07-13 | 39 |
| RVB | 2016 | During | 2016-07-14 | 2016-08-01 | 4 |
| RVB | 2016 | After | 2016-08-02 | 2016-09-16 | 37 |
| RVB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RVB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RVB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RVB | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RVB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RVB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RVB | 2019 | Before | 2019-07-11 | 2019-08-25 | 39 |
| RVB | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RVB | 2019 | After | 2019-09-22 | 2019-11-06 | 36 |
Models
GAM Model 1: Categorical predictors
Plots
df_chla_wt_q %>%
ggplot(aes(x = FlowActionPeriod, y = Chla, fill = FlowActionPeriod)) +
geom_boxplot() +
facet_grid(rows = vars(Year), cols = vars(StationCode)) +
scale_y_log10(labels = trans_format("log10", math_format(10^.x))) +
annotation_logticks(sides = "l") +
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "top"
)
Initial Model
We’ll try running a GAM including all two-way interactions of the categorical predictors - Year, Flow Action Period, Station - and a smooth term for day of year to account for seasonality (restricting the k-value to 5 to reduce overfitting). First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_cat <- gam(
Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)Lets look at the model summary and diagnostics:
summary(m_gam_cat)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY,
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.92533 0.08173 109.206 < 2e-16 ***
## Year_fct2016 -0.12304 0.09641 -1.276 0.201987
## Year_fct2017 0.43975 0.08230 5.343 9.72e-08 ***
## Year_fct2018 0.15352 0.08074 1.901 0.057351 .
## Year_fct2019 0.06346 0.08131 0.781 0.435138
## FlowActionPeriodDuring -0.01236 0.06708 -0.184 0.853795
## FlowActionPeriodAfter 0.58644 0.07712 7.605 3.67e-14 ***
## StationCodeRD22 0.70108 0.08498 8.250 2.23e-16 ***
## StationCodeI80 0.84468 0.08759 9.644 < 2e-16 ***
## StationCodeLIS -0.43029 0.08360 -5.147 2.80e-07 ***
## StationCodeSTTD -0.68280 0.08732 -7.820 7.00e-15 ***
## StationCodeLIB -1.25628 0.08535 -14.720 < 2e-16 ***
## StationCodeRVB -1.28890 0.08343 -15.449 < 2e-16 ***
## Year_fct2016:FlowActionPeriodDuring -0.11791 0.07122 -1.656 0.097897 .
## Year_fct2017:FlowActionPeriodDuring -0.07107 0.05402 -1.316 0.188334
## Year_fct2018:FlowActionPeriodDuring -0.03788 0.05025 -0.754 0.451023
## Year_fct2019:FlowActionPeriodDuring -0.48857 0.05151 -9.485 < 2e-16 ***
## Year_fct2016:FlowActionPeriodAfter -0.16200 0.07631 -2.123 0.033843 *
## Year_fct2017:FlowActionPeriodAfter -0.26360 0.04869 -5.414 6.59e-08 ***
## Year_fct2018:FlowActionPeriodAfter -0.15729 0.04867 -3.232 0.001241 **
## Year_fct2019:FlowActionPeriodAfter 0.06520 0.04933 1.322 0.186312
## Year_fct2016:StationCodeRD22 -0.27091 0.09536 -2.841 0.004527 **
## Year_fct2017:StationCodeRD22 -0.60584 0.09026 -6.713 2.23e-11 ***
## Year_fct2018:StationCodeRD22 -0.42354 0.08782 -4.823 1.48e-06 ***
## Year_fct2019:StationCodeRD22 -0.21418 0.08845 -2.421 0.015512 *
## Year_fct2016:StationCodeI80 0.35595 0.09680 3.677 0.000240 ***
## Year_fct2017:StationCodeI80 -0.23252 0.09193 -2.529 0.011478 *
## Year_fct2018:StationCodeI80 -0.49323 0.08953 -5.509 3.88e-08 ***
## Year_fct2019:StationCodeI80 0.45168 0.09016 5.010 5.73e-07 ***
## Year_fct2016:StationCodeLIS 0.75033 0.09371 8.007 1.60e-15 ***
## Year_fct2017:StationCodeLIS 0.23373 0.09054 2.581 0.009882 **
## Year_fct2018:StationCodeLIS -0.08130 0.08738 -0.930 0.352235
## Year_fct2019:StationCodeLIS 0.63379 0.08855 7.157 1.00e-12 ***
## Year_fct2016:StationCodeSTTD 0.68029 0.09599 7.087 1.66e-12 ***
## Year_fct2017:StationCodeSTTD -0.06767 0.10052 -0.673 0.500848
## Year_fct2018:StationCodeSTTD -0.61835 0.09139 -6.766 1.55e-11 ***
## Year_fct2019:StationCodeSTTD -0.82420 0.08995 -9.163 < 2e-16 ***
## Year_fct2016:StationCodeLIB 1.17255 0.09423 12.443 < 2e-16 ***
## Year_fct2017:StationCodeLIB -0.27805 0.09103 -3.054 0.002272 **
## Year_fct2018:StationCodeLIB -1.92894 0.10180 -18.949 < 2e-16 ***
## Year_fct2019:StationCodeLIB -0.61755 0.10553 -5.852 5.32e-09 ***
## Year_fct2016:StationCodeRVB 0.35867 0.09694 3.700 0.000219 ***
## Year_fct2017:StationCodeRVB -0.55956 0.08929 -6.267 4.14e-10 ***
## Year_fct2018:StationCodeRVB -0.38660 0.08685 -4.452 8.80e-06 ***
## Year_fct2019:StationCodeRVB -0.65441 0.08842 -7.401 1.70e-13 ***
## FlowActionPeriodDuring:StationCodeRD22 -0.41717 0.06339 -6.581 5.40e-11 ***
## FlowActionPeriodAfter:StationCodeRD22 -0.62430 0.05636 -11.077 < 2e-16 ***
## FlowActionPeriodDuring:StationCodeI80 -0.51777 0.06421 -8.063 1.02e-15 ***
## FlowActionPeriodAfter:StationCodeI80 -0.50146 0.05684 -8.822 < 2e-16 ***
## FlowActionPeriodDuring:StationCodeLIS 0.42332 0.06385 6.630 3.90e-11 ***
## FlowActionPeriodAfter:StationCodeLIS -0.52929 0.05589 -9.470 < 2e-16 ***
## FlowActionPeriodDuring:StationCodeSTTD 1.12064 0.06775 16.541 < 2e-16 ***
## FlowActionPeriodAfter:StationCodeSTTD -0.38714 0.06322 -6.123 1.02e-09 ***
## FlowActionPeriodDuring:StationCodeLIB 0.26405 0.06912 3.820 0.000136 ***
## FlowActionPeriodAfter:StationCodeLIB -0.75378 0.06336 -11.897 < 2e-16 ***
## FlowActionPeriodDuring:StationCodeRVB 0.08567 0.06416 1.335 0.181922
## FlowActionPeriodAfter:StationCodeRVB -0.69367 0.05617 -12.350 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.772 3.974 27.07 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.864 Deviance explained = 86.6%
## -REML = 1561.7 Scale est. = 0.13479 n = 3478appraise(m_gam_cat)
k.check(m_gam_cat)## k' edf k-index p-value
## s(DOY) 4 3.77179 0.9457395 0draw(m_gam_cat, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_cat, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_cat))
Box.test(residuals(m_gam_cat), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_cat)
## X-squared = 8609.3, df = 20, p-value < 2.2e-16Model with AR() correlation structure
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_cat <- as.formula("Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, k = 5)")
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_cat <- gamm(
f_gam_cat,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_cat_ar1 <- gamm(
f_gam_cat,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat_ar2 <- gamm(
f_gam_cat,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat_ar3 <- gamm(
f_gam_cat,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_cat$lme,
m_gamm_cat_ar1$lme,
m_gamm_cat_ar2$lme,
m_gamm_cat_ar3$lme
)## Model df AIC BIC logLik Test L.Ratio
## m_gamm_cat$lme 1 60 3243.372 3611.616 -1561.686
## m_gamm_cat_ar1$lme 2 61 -2167.534 -1793.153 1144.767 1 vs 2 5412.906
## m_gamm_cat_ar2$lme 3 62 -2189.412 -1808.893 1156.706 2 vs 3 23.878
## m_gamm_cat_ar3$lme 4 63 -2188.415 -1801.759 1157.208 3 vs 4 1.004
## p-value
## m_gamm_cat$lme
## m_gamm_cat_ar1$lme <.0001
## m_gamm_cat_ar2$lme <.0001
## m_gamm_cat_ar3$lme 0.3164It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
AR(2) Model
# Pull out the residuals and the GAM model
resid_cat_ar2 <- residuals(m_gamm_cat_ar2$lme, type = "normalized")
m_gamm_cat_ar2_gam <- m_gamm_cat_ar2$gam
summary(m_gamm_cat_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY,
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.092494 0.347301 26.180 < 2e-16
## Year_fct2016 0.074100 0.454595 0.163 0.870526
## Year_fct2017 0.351443 0.428540 0.820 0.412220
## Year_fct2018 0.220090 0.419794 0.524 0.600117
## Year_fct2019 -0.004065 0.421429 -0.010 0.992305
## FlowActionPeriodDuring 0.106229 0.094542 1.124 0.261253
## FlowActionPeriodAfter 0.205035 0.143274 1.431 0.152500
## StationCodeRD22 0.213335 0.427422 0.499 0.617727
## StationCodeI80 0.541690 0.439999 1.231 0.218364
## StationCodeLIS -0.597751 0.419022 -1.427 0.153805
## StationCodeSTTD -0.309144 0.429632 -0.720 0.471848
## StationCodeLIB -1.392380 0.421407 -3.304 0.000963
## StationCodeRVB -1.450049 0.417001 -3.477 0.000513
## Year_fct2016:FlowActionPeriodDuring -0.162149 0.090727 -1.787 0.073991
## Year_fct2017:FlowActionPeriodDuring 0.011836 0.089175 0.133 0.894417
## Year_fct2018:FlowActionPeriodDuring -0.171794 0.089080 -1.929 0.053872
## Year_fct2019:FlowActionPeriodDuring 0.054852 0.089154 0.615 0.538428
## Year_fct2016:FlowActionPeriodAfter -0.216376 0.127679 -1.695 0.090227
## Year_fct2017:FlowActionPeriodAfter -0.082755 0.125185 -0.661 0.508619
## Year_fct2018:FlowActionPeriodAfter -0.235315 0.125043 -1.882 0.059938
## Year_fct2019:FlowActionPeriodAfter 0.127736 0.124974 1.022 0.306806
## Year_fct2016:StationCodeRD22 -0.409097 0.576632 -0.709 0.478087
## Year_fct2017:StationCodeRD22 -0.510460 0.546951 -0.933 0.350740
## Year_fct2018:StationCodeRD22 -0.375917 0.536410 -0.701 0.483474
## Year_fct2019:StationCodeRD22 -0.125667 0.538556 -0.233 0.815511
## Year_fct2016:StationCodeI80 0.159128 0.585194 0.272 0.785697
## Year_fct2017:StationCodeI80 -0.229520 0.556103 -0.413 0.679830
## Year_fct2018:StationCodeI80 -0.577356 0.545768 -1.058 0.290186
## Year_fct2019:StationCodeI80 0.412105 0.547880 0.752 0.451994
## Year_fct2016:StationCodeLIS 0.766734 0.561958 1.364 0.172532
## Year_fct2017:StationCodeLIS 0.178449 0.544443 0.328 0.743110
## Year_fct2018:StationCodeLIS -0.130811 0.531634 -0.246 0.805655
## Year_fct2019:StationCodeLIS 0.384296 0.537308 0.715 0.474519
## Year_fct2016:StationCodeSTTD 0.494794 0.579473 0.854 0.393237
## Year_fct2017:StationCodeSTTD -0.666965 0.584959 -1.140 0.254287
## Year_fct2018:StationCodeSTTD -0.891467 0.556482 -1.602 0.109255
## Year_fct2019:StationCodeSTTD -1.068222 0.547691 -1.950 0.051209
## Year_fct2016:StationCodeLIB 1.032582 0.565695 1.825 0.068038
## Year_fct2017:StationCodeLIB -0.298959 0.549421 -0.544 0.586385
## Year_fct2018:StationCodeLIB -2.012752 0.597348 -3.369 0.000761
## Year_fct2019:StationCodeLIB -0.761233 0.605732 -1.257 0.208942
## Year_fct2016:StationCodeRVB 0.385100 0.575289 0.669 0.503284
## Year_fct2017:StationCodeRVB -0.578913 0.539279 -1.073 0.283125
## Year_fct2018:StationCodeRVB -0.447074 0.528583 -0.846 0.397726
## Year_fct2019:StationCodeRVB -0.854363 0.537540 -1.589 0.112064
## FlowActionPeriodDuring:StationCodeRD22 -0.067302 0.105406 -0.638 0.523192
## FlowActionPeriodAfter:StationCodeRD22 0.031609 0.148699 0.213 0.831675
## FlowActionPeriodDuring:StationCodeI80 -0.058480 0.105829 -0.553 0.580579
## FlowActionPeriodAfter:StationCodeI80 -0.002189 0.149048 -0.015 0.988284
## FlowActionPeriodDuring:StationCodeLIS -0.087621 0.105289 -0.832 0.405357
## FlowActionPeriodAfter:StationCodeLIS 0.013527 0.148521 0.091 0.927435
## FlowActionPeriodDuring:StationCodeSTTD 0.056894 0.105964 0.537 0.591360
## FlowActionPeriodAfter:StationCodeSTTD -0.493124 0.150394 -3.279 0.001053
## FlowActionPeriodDuring:StationCodeLIB -0.098209 0.106665 -0.921 0.357264
## FlowActionPeriodAfter:StationCodeLIB -0.224653 0.150868 -1.489 0.136561
## FlowActionPeriodDuring:StationCodeRVB -0.067906 0.105332 -0.645 0.519175
## FlowActionPeriodAfter:StationCodeRVB -0.174298 0.148565 -1.173 0.240792
##
## (Intercept) ***
## Year_fct2016
## Year_fct2017
## Year_fct2018
## Year_fct2019
## FlowActionPeriodDuring
## FlowActionPeriodAfter
## StationCodeRD22
## StationCodeI80
## StationCodeLIS
## StationCodeSTTD
## StationCodeLIB ***
## StationCodeRVB ***
## Year_fct2016:FlowActionPeriodDuring .
## Year_fct2017:FlowActionPeriodDuring
## Year_fct2018:FlowActionPeriodDuring .
## Year_fct2019:FlowActionPeriodDuring
## Year_fct2016:FlowActionPeriodAfter .
## Year_fct2017:FlowActionPeriodAfter
## Year_fct2018:FlowActionPeriodAfter .
## Year_fct2019:FlowActionPeriodAfter
## Year_fct2016:StationCodeRD22
## Year_fct2017:StationCodeRD22
## Year_fct2018:StationCodeRD22
## Year_fct2019:StationCodeRD22
## Year_fct2016:StationCodeI80
## Year_fct2017:StationCodeI80
## Year_fct2018:StationCodeI80
## Year_fct2019:StationCodeI80
## Year_fct2016:StationCodeLIS
## Year_fct2017:StationCodeLIS
## Year_fct2018:StationCodeLIS
## Year_fct2019:StationCodeLIS
## Year_fct2016:StationCodeSTTD
## Year_fct2017:StationCodeSTTD
## Year_fct2018:StationCodeSTTD
## Year_fct2019:StationCodeSTTD .
## Year_fct2016:StationCodeLIB .
## Year_fct2017:StationCodeLIB
## Year_fct2018:StationCodeLIB ***
## Year_fct2019:StationCodeLIB
## Year_fct2016:StationCodeRVB
## Year_fct2017:StationCodeRVB
## Year_fct2018:StationCodeRVB
## Year_fct2019:StationCodeRVB
## FlowActionPeriodDuring:StationCodeRD22
## FlowActionPeriodAfter:StationCodeRD22
## FlowActionPeriodDuring:StationCodeI80
## FlowActionPeriodAfter:StationCodeI80
## FlowActionPeriodDuring:StationCodeLIS
## FlowActionPeriodAfter:StationCodeLIS
## FlowActionPeriodDuring:StationCodeSTTD
## FlowActionPeriodAfter:StationCodeSTTD **
## FlowActionPeriodDuring:StationCodeLIB
## FlowActionPeriodAfter:StationCodeLIB
## FlowActionPeriodDuring:StationCodeRVB
## FlowActionPeriodAfter:StationCodeRVB
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.734 3.734 11.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.785
## Scale est. = 0.27538 n = 3478appraise(m_gamm_cat_ar2_gam)
k.check(m_gamm_cat_ar2_gam)## k' edf k-index p-value
## s(DOY) 4 3.734167 0.7698425 0draw(m_gamm_cat_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_cat_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_cat_ar2)
Box.test(resid_cat_ar2, lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: resid_cat_ar2
## X-squared = 32.886, df = 20, p-value = 0.03473The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_cat_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY,
## k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 0.345 0.8477
## FlowActionPeriod 2 1.070 0.3432
## StationCode 6 9.384 3.08e-10
## Year_fct:FlowActionPeriod 8 1.991 0.0438
## Year_fct:StationCode 24 2.732 1.16e-05
## FlowActionPeriod:StationCode 12 4.817 6.68e-08
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.734 3.734 11.49 <2e-16All the interaction terms are significant in this model, so we’ll use this one for the model selection process.
GAM Model 2: Categorical predictors adding water temperature
Initial Model
We’ll try running the categorical model adding water temperature to see if it improves the predictive power of the model.
m_gam_cat_wt <- gam(
Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)Lets look at the model summary and diagnostics:
summary(m_gam_cat_wt)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 +
## s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.4732371 0.4077027 23.236 < 2e-16
## Year_fct2016 0.1185304 0.3872590 0.306 0.759566
## Year_fct2017 0.1743035 0.3365746 0.518 0.604579
## Year_fct2018 0.3038120 0.2848691 1.066 0.286275
## Year_fct2019 0.0558398 0.2655069 0.210 0.833435
## FlowActionPeriodDuring -0.8448474 0.3953310 -2.137 0.032664
## FlowActionPeriodAfter 0.1409841 0.3375317 0.418 0.676199
## StationCodeRD22 -0.0797142 0.3149937 -0.253 0.800233
## StationCodeI80 1.0484912 0.3010805 3.482 0.000503
## StationCodeLIS -3.0959341 0.2998499 -10.325 < 2e-16
## StationCodeSTTD -2.1712782 0.3221696 -6.740 1.86e-11
## StationCodeLIB -2.4195633 0.3868708 -6.254 4.49e-10
## StationCodeRVB -2.5192484 0.3578309 -7.040 2.31e-12
## WaterTemp -0.0172525 0.0159383 -1.082 0.279127
## Year_fct2016:FlowActionPeriodDuring -0.1478412 0.0717935 -2.059 0.039545
## Year_fct2017:FlowActionPeriodDuring -0.1142243 0.0559969 -2.040 0.041444
## Year_fct2018:FlowActionPeriodDuring 0.0081453 0.0566778 0.144 0.885736
## Year_fct2019:FlowActionPeriodDuring -0.4834680 0.0537729 -8.991 < 2e-16
## Year_fct2016:FlowActionPeriodAfter -0.0463716 0.0939445 -0.494 0.621616
## Year_fct2017:FlowActionPeriodAfter -0.1164966 0.0985549 -1.182 0.237269
## Year_fct2018:FlowActionPeriodAfter -0.1779551 0.0818456 -2.174 0.029753
## Year_fct2019:FlowActionPeriodAfter 0.1434797 0.0868389 1.652 0.098575
## Year_fct2016:StationCodeRD22 -0.4030850 0.1061711 -3.797 0.000149
## Year_fct2017:StationCodeRD22 -0.6258445 0.0886093 -7.063 1.97e-12
## Year_fct2018:StationCodeRD22 -0.4032341 0.0861352 -4.681 2.96e-06
## Year_fct2019:StationCodeRD22 -0.2105029 0.0863540 -2.438 0.014833
## Year_fct2016:StationCodeI80 0.3764403 0.1078367 3.491 0.000488
## Year_fct2017:StationCodeI80 -0.2152108 0.0907540 -2.371 0.017778
## Year_fct2018:StationCodeI80 -0.5036823 0.0880693 -5.719 1.16e-08
## Year_fct2019:StationCodeI80 0.4707630 0.0885154 5.318 1.11e-07
## Year_fct2016:StationCodeLIS 0.3980072 0.1037768 3.835 0.000128
## Year_fct2017:StationCodeLIS 0.1331875 0.0890675 1.495 0.134915
## Year_fct2018:StationCodeLIS -0.0665137 0.0856235 -0.777 0.437321
## Year_fct2019:StationCodeLIS 0.5866469 0.0867510 6.762 1.59e-11
## Year_fct2016:StationCodeSTTD 0.4647981 0.1075972 4.320 1.61e-05
## Year_fct2017:StationCodeSTTD -0.1040733 0.0993597 -1.047 0.294971
## Year_fct2018:StationCodeSTTD -0.5941584 0.0894294 -6.644 3.54e-11
## Year_fct2019:StationCodeSTTD -0.8406972 0.0880556 -9.547 < 2e-16
## Year_fct2016:StationCodeLIB 1.0594550 0.1104642 9.591 < 2e-16
## Year_fct2017:StationCodeLIB -0.2288879 0.0930776 -2.459 0.013978
## Year_fct2018:StationCodeLIB -1.8694287 0.1011190 -18.487 < 2e-16
## Year_fct2019:StationCodeLIB -0.5687776 0.1049190 -5.421 6.34e-08
## Year_fct2016:StationCodeRVB 0.2273758 0.1104979 2.058 0.039691
## Year_fct2017:StationCodeRVB -0.4993343 0.0911536 -5.478 4.62e-08
## Year_fct2018:StationCodeRVB -0.3330771 0.0868161 -3.837 0.000127
## Year_fct2019:StationCodeRVB -0.5914447 0.0887079 -6.667 3.03e-11
## Year_fct2016:WaterTemp -0.0047186 0.0147895 -0.319 0.749709
## Year_fct2017:WaterTemp 0.0102785 0.0131868 0.779 0.435767
## Year_fct2018:WaterTemp -0.0074263 0.0113164 -0.656 0.511709
## Year_fct2019:WaterTemp -0.0009458 0.0102805 -0.092 0.926706
## FlowActionPeriodDuring:StationCodeRD22 -0.3367835 0.0692219 -4.865 1.19e-06
## FlowActionPeriodAfter:StationCodeRD22 -0.3952011 0.1035785 -3.815 0.000138
## FlowActionPeriodDuring:StationCodeI80 -0.5058395 0.0700402 -7.222 6.28e-13
## FlowActionPeriodAfter:StationCodeI80 -0.5403872 0.1029498 -5.249 1.62e-07
## FlowActionPeriodDuring:StationCodeLIS 0.6066616 0.0696598 8.709 < 2e-16
## FlowActionPeriodAfter:StationCodeLIS 0.1670770 0.0997364 1.675 0.093989
## FlowActionPeriodDuring:StationCodeSTTD 1.2258165 0.0734088 16.699 < 2e-16
## FlowActionPeriodAfter:StationCodeSTTD 0.0006864 0.1043389 0.007 0.994751
## FlowActionPeriodDuring:StationCodeLIB 0.4149955 0.0840477 4.938 8.29e-07
## FlowActionPeriodAfter:StationCodeLIB -0.4732155 0.1097713 -4.311 1.67e-05
## FlowActionPeriodDuring:StationCodeRVB 0.2264187 0.0775096 2.921 0.003510
## FlowActionPeriodAfter:StationCodeRVB -0.3932928 0.1064178 -3.696 0.000223
## FlowActionPeriodDuring:WaterTemp 0.0304368 0.0159083 1.913 0.055797
## FlowActionPeriodAfter:WaterTemp 0.0026541 0.0125480 0.212 0.832498
## StationCodeRD22:WaterTemp 0.0307775 0.0121747 2.528 0.011517
## StationCodeI80:WaterTemp -0.0098427 0.0115431 -0.853 0.393890
## StationCodeLIS:WaterTemp 0.1106623 0.0116009 9.539 < 2e-16
## StationCodeSTTD:WaterTemp 0.0615658 0.0127363 4.834 1.40e-06
## StationCodeLIB:WaterTemp 0.0474936 0.0159554 2.977 0.002935
## StationCodeRVB:WaterTemp 0.0499363 0.0143724 3.474 0.000518
##
## (Intercept) ***
## Year_fct2016
## Year_fct2017
## Year_fct2018
## Year_fct2019
## FlowActionPeriodDuring *
## FlowActionPeriodAfter
## StationCodeRD22
## StationCodeI80 ***
## StationCodeLIS ***
## StationCodeSTTD ***
## StationCodeLIB ***
## StationCodeRVB ***
## WaterTemp
## Year_fct2016:FlowActionPeriodDuring *
## Year_fct2017:FlowActionPeriodDuring *
## Year_fct2018:FlowActionPeriodDuring
## Year_fct2019:FlowActionPeriodDuring ***
## Year_fct2016:FlowActionPeriodAfter
## Year_fct2017:FlowActionPeriodAfter
## Year_fct2018:FlowActionPeriodAfter *
## Year_fct2019:FlowActionPeriodAfter .
## Year_fct2016:StationCodeRD22 ***
## Year_fct2017:StationCodeRD22 ***
## Year_fct2018:StationCodeRD22 ***
## Year_fct2019:StationCodeRD22 *
## Year_fct2016:StationCodeI80 ***
## Year_fct2017:StationCodeI80 *
## Year_fct2018:StationCodeI80 ***
## Year_fct2019:StationCodeI80 ***
## Year_fct2016:StationCodeLIS ***
## Year_fct2017:StationCodeLIS
## Year_fct2018:StationCodeLIS
## Year_fct2019:StationCodeLIS ***
## Year_fct2016:StationCodeSTTD ***
## Year_fct2017:StationCodeSTTD
## Year_fct2018:StationCodeSTTD ***
## Year_fct2019:StationCodeSTTD ***
## Year_fct2016:StationCodeLIB ***
## Year_fct2017:StationCodeLIB *
## Year_fct2018:StationCodeLIB ***
## Year_fct2019:StationCodeLIB ***
## Year_fct2016:StationCodeRVB *
## Year_fct2017:StationCodeRVB ***
## Year_fct2018:StationCodeRVB ***
## Year_fct2019:StationCodeRVB ***
## Year_fct2016:WaterTemp
## Year_fct2017:WaterTemp
## Year_fct2018:WaterTemp
## Year_fct2019:WaterTemp
## FlowActionPeriodDuring:StationCodeRD22 ***
## FlowActionPeriodAfter:StationCodeRD22 ***
## FlowActionPeriodDuring:StationCodeI80 ***
## FlowActionPeriodAfter:StationCodeI80 ***
## FlowActionPeriodDuring:StationCodeLIS ***
## FlowActionPeriodAfter:StationCodeLIS .
## FlowActionPeriodDuring:StationCodeSTTD ***
## FlowActionPeriodAfter:StationCodeSTTD
## FlowActionPeriodDuring:StationCodeLIB ***
## FlowActionPeriodAfter:StationCodeLIB ***
## FlowActionPeriodDuring:StationCodeRVB **
## FlowActionPeriodAfter:StationCodeRVB ***
## FlowActionPeriodDuring:WaterTemp .
## FlowActionPeriodAfter:WaterTemp
## StationCodeRD22:WaterTemp *
## StationCodeI80:WaterTemp
## StationCodeLIS:WaterTemp ***
## StationCodeSTTD:WaterTemp ***
## StationCodeLIB:WaterTemp **
## StationCodeRVB:WaterTemp ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.731 3.966 23.94 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.87 Deviance explained = 87.3%
## -REML = 1517.2 Scale est. = 0.12822 n = 3478appraise(m_gam_cat_wt)
k.check(m_gam_cat_wt)## k' edf k-index p-value
## s(DOY) 4 3.730833 0.9463198 0draw(m_gam_cat_wt, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_cat_wt, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_cat_wt))
Box.test(residuals(m_gam_cat_wt), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_cat_wt)
## X-squared = 8031.8, df = 20, p-value < 2.2e-16Model with AR() correlation structure
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_cat_wt <- as.formula("Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 + s(DOY, k = 5)")
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_cat_wt <- gamm(
f_gam_cat_wt,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_cat_wt_ar1 <- gamm(
f_gam_cat_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat_wt_ar2 <- gamm(
f_gam_cat_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat_wt_ar3 <- gamm(
f_gam_cat_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_cat_wt$lme,
m_gamm_cat_wt_ar1$lme,
m_gamm_cat_wt_ar2$lme,
m_gamm_cat_wt_ar3$lme
)## Model df AIC BIC logLik Test L.Ratio
## m_gamm_cat_wt$lme 1 73 3180.470 3628.222 -1517.235
## m_gamm_cat_wt_ar1$lme 2 74 -2135.571 -1681.686 1141.786 1 vs 2 5318.041
## m_gamm_cat_wt_ar2$lme 3 75 -2160.486 -1700.467 1155.243 2 vs 3 26.915
## m_gamm_cat_wt_ar3$lme 4 76 -2159.702 -1693.550 1155.851 3 vs 4 1.216
## p-value
## m_gamm_cat_wt$lme
## m_gamm_cat_wt_ar1$lme <.0001
## m_gamm_cat_wt_ar2$lme <.0001
## m_gamm_cat_wt_ar3$lme 0.2701It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
AR(2) Model
# Pull out the residuals and the GAM model
resid_cat_wt_ar2 <- residuals(m_gamm_cat_wt_ar2$lme, type = "normalized")
m_gamm_cat_wt_ar2_gam <- m_gamm_cat_wt_ar2$gam
summary(m_gamm_cat_wt_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 +
## s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.1445090 0.6901525 13.250 < 2e-16
## Year_fct2016 1.3780976 0.6921213 1.991 0.04655
## Year_fct2017 0.0579152 0.6440942 0.090 0.92836
## Year_fct2018 -0.3370015 0.6545952 -0.515 0.60671
## Year_fct2019 -0.3900285 0.6187327 -0.630 0.52850
## FlowActionPeriodDuring -0.4663337 0.4240927 -1.100 0.27158
## FlowActionPeriodAfter -0.3836823 0.3632026 -1.056 0.29087
## StationCodeRD22 -1.1868159 0.6842598 -1.734 0.08293
## StationCodeI80 0.9164038 0.6678780 1.372 0.17012
## StationCodeLIS -1.6220417 0.6764787 -2.398 0.01655
## StationCodeSTTD -0.5586425 0.6987943 -0.799 0.42409
## StationCodeLIB -0.2983707 0.7674643 -0.389 0.69747
## StationCodeRVB -1.9098221 0.8359931 -2.284 0.02240
## WaterTemp 0.0009966 0.0243656 0.041 0.96738
## Year_fct2016:FlowActionPeriodDuring -0.1572496 0.0926170 -1.698 0.08963
## Year_fct2017:FlowActionPeriodDuring -0.0258551 0.0919774 -0.281 0.77865
## Year_fct2018:FlowActionPeriodDuring -0.1036021 0.0974311 -1.063 0.28770
## Year_fct2019:FlowActionPeriodDuring 0.0270928 0.0904924 0.299 0.76466
## Year_fct2016:FlowActionPeriodAfter -0.1711207 0.1383100 -1.237 0.21609
## Year_fct2017:FlowActionPeriodAfter -0.0309569 0.1287914 -0.240 0.81006
## Year_fct2018:FlowActionPeriodAfter -0.1354965 0.1343868 -1.008 0.31340
## Year_fct2019:FlowActionPeriodAfter 0.1469567 0.1263950 1.163 0.24504
## Year_fct2016:StationCodeRD22 -0.6182502 0.6279173 -0.985 0.32489
## Year_fct2017:StationCodeRD22 -0.5837985 0.5953210 -0.981 0.32684
## Year_fct2018:StationCodeRD22 -0.3239681 0.5851972 -0.554 0.57989
## Year_fct2019:StationCodeRD22 -0.1547261 0.5869605 -0.264 0.79210
## Year_fct2016:StationCodeI80 0.2258097 0.6387648 0.354 0.72373
## Year_fct2017:StationCodeI80 -0.1772145 0.6064542 -0.292 0.77014
## Year_fct2018:StationCodeI80 -0.5434617 0.5963020 -0.911 0.36216
## Year_fct2019:StationCodeI80 0.4333058 0.5981682 0.724 0.46888
## Year_fct2016:StationCodeLIS 0.5871339 0.6122516 0.959 0.33764
## Year_fct2017:StationCodeLIS 0.1547944 0.5925727 0.261 0.79394
## Year_fct2018:StationCodeLIS -0.1245806 0.5799582 -0.215 0.82993
## Year_fct2019:StationCodeLIS 0.3524021 0.5857178 0.602 0.54744
## Year_fct2016:StationCodeSTTD 0.4354888 0.6320048 0.689 0.49083
## Year_fct2017:StationCodeSTTD -0.6528634 0.6368032 -1.025 0.30533
## Year_fct2018:StationCodeSTTD -0.8898961 0.6067842 -1.467 0.14258
## Year_fct2019:StationCodeSTTD -1.0781623 0.5972945 -1.805 0.07115
## Year_fct2016:StationCodeLIB 0.9323841 0.6169532 1.511 0.13081
## Year_fct2017:StationCodeLIB -0.2392731 0.5986330 -0.400 0.68940
## Year_fct2018:StationCodeLIB -2.0226118 0.6499257 -3.112 0.00187
## Year_fct2019:StationCodeLIB -0.7318547 0.6572217 -1.114 0.26555
## Year_fct2016:StationCodeRVB 0.2238684 0.6270456 0.357 0.72110
## Year_fct2017:StationCodeRVB -0.5138895 0.5877741 -0.874 0.38202
## Year_fct2018:StationCodeRVB -0.3928585 0.5778469 -0.680 0.49664
## Year_fct2019:StationCodeRVB -0.7921245 0.5874493 -1.348 0.17762
## Year_fct2016:WaterTemp -0.0524611 0.0207223 -2.532 0.01140
## Year_fct2017:WaterTemp 0.0126743 0.0191452 0.662 0.50801
## Year_fct2018:WaterTemp 0.0245102 0.0205942 1.190 0.23407
## Year_fct2019:WaterTemp 0.0193582 0.0181899 1.064 0.28730
## FlowActionPeriodDuring:StationCodeRD22 -0.0578678 0.1049823 -0.551 0.58152
## FlowActionPeriodAfter:StationCodeRD22 0.0980816 0.1507379 0.651 0.51530
## FlowActionPeriodDuring:StationCodeI80 -0.0674505 0.1057957 -0.638 0.52381
## FlowActionPeriodAfter:StationCodeI80 -0.0037150 0.1519709 -0.024 0.98050
## FlowActionPeriodDuring:StationCodeLIS -0.0502238 0.1055792 -0.476 0.63432
## FlowActionPeriodAfter:StationCodeLIS 0.1157942 0.1515187 0.764 0.44479
## FlowActionPeriodDuring:StationCodeSTTD 0.0666638 0.1064582 0.626 0.53123
## FlowActionPeriodAfter:StationCodeSTTD -0.4517271 0.1532547 -2.948 0.00322
## FlowActionPeriodDuring:StationCodeLIB -0.0730926 0.1151386 -0.635 0.52559
## FlowActionPeriodAfter:StationCodeLIB -0.2143517 0.1567062 -1.368 0.17145
## FlowActionPeriodDuring:StationCodeRVB -0.0110737 0.1127300 -0.098 0.92175
## FlowActionPeriodAfter:StationCodeRVB -0.0964396 0.1546565 -0.624 0.53295
## FlowActionPeriodDuring:WaterTemp 0.0231542 0.0168338 1.375 0.16908
## FlowActionPeriodAfter:WaterTemp 0.0210524 0.0143386 1.468 0.14213
## StationCodeRD22:WaterTemp 0.0653478 0.0222444 2.938 0.00333
## StationCodeI80:WaterTemp -0.0191793 0.0207291 -0.925 0.35491
## StationCodeLIS:WaterTemp 0.0464759 0.0218848 2.124 0.03377
## StationCodeSTTD:WaterTemp 0.0103409 0.0232382 0.445 0.65635
## StationCodeLIB:WaterTemp -0.0551245 0.0276723 -1.992 0.04645
## StationCodeRVB:WaterTemp 0.0192565 0.0319521 0.603 0.54677
##
## (Intercept) ***
## Year_fct2016 *
## Year_fct2017
## Year_fct2018
## Year_fct2019
## FlowActionPeriodDuring
## FlowActionPeriodAfter
## StationCodeRD22 .
## StationCodeI80
## StationCodeLIS *
## StationCodeSTTD
## StationCodeLIB
## StationCodeRVB *
## WaterTemp
## Year_fct2016:FlowActionPeriodDuring .
## Year_fct2017:FlowActionPeriodDuring
## Year_fct2018:FlowActionPeriodDuring
## Year_fct2019:FlowActionPeriodDuring
## Year_fct2016:FlowActionPeriodAfter
## Year_fct2017:FlowActionPeriodAfter
## Year_fct2018:FlowActionPeriodAfter
## Year_fct2019:FlowActionPeriodAfter
## Year_fct2016:StationCodeRD22
## Year_fct2017:StationCodeRD22
## Year_fct2018:StationCodeRD22
## Year_fct2019:StationCodeRD22
## Year_fct2016:StationCodeI80
## Year_fct2017:StationCodeI80
## Year_fct2018:StationCodeI80
## Year_fct2019:StationCodeI80
## Year_fct2016:StationCodeLIS
## Year_fct2017:StationCodeLIS
## Year_fct2018:StationCodeLIS
## Year_fct2019:StationCodeLIS
## Year_fct2016:StationCodeSTTD
## Year_fct2017:StationCodeSTTD
## Year_fct2018:StationCodeSTTD
## Year_fct2019:StationCodeSTTD .
## Year_fct2016:StationCodeLIB
## Year_fct2017:StationCodeLIB
## Year_fct2018:StationCodeLIB **
## Year_fct2019:StationCodeLIB
## Year_fct2016:StationCodeRVB
## Year_fct2017:StationCodeRVB
## Year_fct2018:StationCodeRVB
## Year_fct2019:StationCodeRVB
## Year_fct2016:WaterTemp *
## Year_fct2017:WaterTemp
## Year_fct2018:WaterTemp
## Year_fct2019:WaterTemp
## FlowActionPeriodDuring:StationCodeRD22
## FlowActionPeriodAfter:StationCodeRD22
## FlowActionPeriodDuring:StationCodeI80
## FlowActionPeriodAfter:StationCodeI80
## FlowActionPeriodDuring:StationCodeLIS
## FlowActionPeriodAfter:StationCodeLIS
## FlowActionPeriodDuring:StationCodeSTTD
## FlowActionPeriodAfter:StationCodeSTTD **
## FlowActionPeriodDuring:StationCodeLIB
## FlowActionPeriodAfter:StationCodeLIB
## FlowActionPeriodDuring:StationCodeRVB
## FlowActionPeriodAfter:StationCodeRVB
## FlowActionPeriodDuring:WaterTemp
## FlowActionPeriodAfter:WaterTemp
## StationCodeRD22:WaterTemp **
## StationCodeI80:WaterTemp
## StationCodeLIS:WaterTemp *
## StationCodeSTTD:WaterTemp
## StationCodeLIB:WaterTemp *
## StationCodeRVB:WaterTemp
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.669 3.669 8.849 5.92e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.774
## Scale est. = 0.30258 n = 3478appraise(m_gamm_cat_wt_ar2_gam)
k.check(m_gamm_cat_wt_ar2_gam)## k' edf k-index p-value
## s(DOY) 4 3.669285 0.7410218 0draw(m_gamm_cat_wt_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_cat_wt_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_cat_wt_ar2)
Box.test(resid_cat_wt_ar2, lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: resid_cat_wt_ar2
## X-squared = 30.859, df = 20, p-value = 0.05708The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_cat_wt_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 +
## s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 2.412 0.047021
## FlowActionPeriod 2 0.740 0.477314
## StationCode 6 5.258 2.11e-05
## WaterTemp 1 0.002 0.967377
## Year_fct:FlowActionPeriod 8 1.159 0.320505
## Year_fct:StationCode 24 2.285 0.000357
## Year_fct:WaterTemp 4 4.446 0.001386
## FlowActionPeriod:StationCode 12 4.941 3.63e-08
## FlowActionPeriod:WaterTemp 2 1.414 0.243334
## StationCode:WaterTemp 6 6.427 9.53e-07
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.669 3.669 8.849 5.92e-06The Year:StationCode, Year:WaterTemp, FlowActionPeriod:StationCode, and StationCode:WaterTemp interactions were significant among the interaction terms of the model. Let’s remove the Year:FlowActionPeriod and FlowActionPeriod:WaterTemp interaction terms and check the model summary and diagnostics.
m_gamm_cat_wt_ar2b <- gamm(
Chla_log ~ (Year_fct + StationCode + WaterTemp)^2 + StationCode * FlowActionPeriod + s(DOY, k = 5),
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
# Pull out the residuals and the GAM model
resid_cat_wt_ar2b <- residuals(m_gamm_cat_wt_ar2b$lme, type = "normalized")
m_gamm_cat_wt_ar2b_gam <- m_gamm_cat_wt_ar2b$gam
summary(m_gamm_cat_wt_ar2b_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + StationCode + WaterTemp)^2 + StationCode *
## FlowActionPeriod + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.782805 0.664490 13.217 < 2e-16
## Year_fct2016 1.341143 0.707272 1.896 0.05802
## Year_fct2017 -0.006022 0.652080 -0.009 0.99263
## Year_fct2018 -0.528205 0.657513 -0.803 0.42184
## Year_fct2019 -0.272579 0.628678 -0.434 0.66462
## StationCodeRD22 -1.163183 0.707102 -1.645 0.10006
## StationCodeI80 0.951126 0.691075 1.376 0.16882
## StationCodeLIS -1.559134 0.698191 -2.233 0.02561
## StationCodeSTTD -0.493857 0.720065 -0.686 0.49285
## StationCodeLIB -0.085665 0.784497 -0.109 0.91305
## StationCodeRVB -1.782352 0.855977 -2.082 0.03739
## WaterTemp 0.016614 0.022904 0.725 0.46828
## FlowActionPeriodDuring 0.052204 0.074987 0.696 0.48637
## FlowActionPeriodAfter 0.099477 0.110745 0.898 0.36912
## Year_fct2016:StationCodeRD22 -0.613100 0.667140 -0.919 0.35816
## Year_fct2017:StationCodeRD22 -0.578281 0.633452 -0.913 0.36136
## Year_fct2018:StationCodeRD22 -0.315304 0.622783 -0.506 0.61269
## Year_fct2019:StationCodeRD22 -0.149630 0.624636 -0.240 0.81070
## Year_fct2016:StationCodeI80 0.223843 0.678341 0.330 0.74143
## Year_fct2017:StationCodeI80 -0.181161 0.644952 -0.281 0.77881
## Year_fct2018:StationCodeI80 -0.547381 0.634232 -0.863 0.38817
## Year_fct2019:StationCodeI80 0.423004 0.636191 0.665 0.50616
## Year_fct2016:StationCodeLIS 0.619479 0.650406 0.952 0.34094
## Year_fct2017:StationCodeLIS 0.170080 0.630490 0.270 0.78736
## Year_fct2018:StationCodeLIS -0.116806 0.617195 -0.189 0.84991
## Year_fct2019:StationCodeLIS 0.349457 0.623223 0.561 0.57502
## Year_fct2016:StationCodeSTTD 0.445998 0.671203 0.664 0.50643
## Year_fct2017:StationCodeSTTD -0.666675 0.676269 -0.986 0.32429
## Year_fct2018:StationCodeSTTD -0.889459 0.645302 -1.378 0.16818
## Year_fct2019:StationCodeSTTD -1.078389 0.635280 -1.698 0.08969
## Year_fct2016:StationCodeLIB 0.940927 0.655358 1.436 0.15117
## Year_fct2017:StationCodeLIB -0.223119 0.636768 -0.350 0.72607
## Year_fct2018:StationCodeLIB -2.019727 0.689265 -2.930 0.00341
## Year_fct2019:StationCodeLIB -0.746492 0.696748 -1.071 0.28407
## Year_fct2016:StationCodeRVB 0.239251 0.665883 0.359 0.71939
## Year_fct2017:StationCodeRVB -0.506117 0.625396 -0.809 0.41841
## Year_fct2018:StationCodeRVB -0.384654 0.614812 -0.626 0.53159
## Year_fct2019:StationCodeRVB -0.807469 0.624907 -1.292 0.19640
## Year_fct2016:WaterTemp -0.056611 0.020507 -2.761 0.00580
## Year_fct2017:WaterTemp 0.013897 0.018782 0.740 0.45940
## Year_fct2018:WaterTemp 0.029403 0.020156 1.459 0.14472
## Year_fct2019:WaterTemp 0.016339 0.017903 0.913 0.36148
## StationCodeRD22:WaterTemp 0.064117 0.022486 2.851 0.00438
## StationCodeI80:WaterTemp -0.020588 0.020915 -0.984 0.32500
## StationCodeLIS:WaterTemp 0.043491 0.022089 1.969 0.04904
## StationCodeSTTD:WaterTemp 0.008117 0.023400 0.347 0.72870
## StationCodeLIB:WaterTemp -0.063750 0.027815 -2.292 0.02197
## StationCodeRVB:WaterTemp 0.015016 0.032392 0.464 0.64298
## StationCodeRD22:FlowActionPeriodDuring -0.054000 0.105338 -0.513 0.60824
## StationCodeI80:FlowActionPeriodDuring -0.062165 0.106067 -0.586 0.55785
## StationCodeLIS:FlowActionPeriodDuring -0.061674 0.105451 -0.585 0.55868
## StationCodeSTTD:FlowActionPeriodDuring 0.042061 0.105724 0.398 0.69078
## StationCodeLIB:FlowActionPeriodDuring -0.133669 0.107077 -1.248 0.21199
## StationCodeRVB:FlowActionPeriodDuring -0.063906 0.105456 -0.606 0.54456
## StationCodeRD22:FlowActionPeriodAfter 0.101187 0.150452 0.673 0.50128
## StationCodeI80:FlowActionPeriodAfter 0.005898 0.150783 0.039 0.96880
## StationCodeLIS:FlowActionPeriodAfter 0.107314 0.151165 0.710 0.47781
## StationCodeSTTD:FlowActionPeriodAfter -0.476380 0.152352 -3.127 0.00178
## StationCodeLIB:FlowActionPeriodAfter -0.265101 0.153142 -1.731 0.08353
## StationCodeRVB:FlowActionPeriodAfter -0.142648 0.151288 -0.943 0.34580
##
## (Intercept) ***
## Year_fct2016 .
## Year_fct2017
## Year_fct2018
## Year_fct2019
## StationCodeRD22
## StationCodeI80
## StationCodeLIS *
## StationCodeSTTD
## StationCodeLIB
## StationCodeRVB *
## WaterTemp
## FlowActionPeriodDuring
## FlowActionPeriodAfter
## Year_fct2016:StationCodeRD22
## Year_fct2017:StationCodeRD22
## Year_fct2018:StationCodeRD22
## Year_fct2019:StationCodeRD22
## Year_fct2016:StationCodeI80
## Year_fct2017:StationCodeI80
## Year_fct2018:StationCodeI80
## Year_fct2019:StationCodeI80
## Year_fct2016:StationCodeLIS
## Year_fct2017:StationCodeLIS
## Year_fct2018:StationCodeLIS
## Year_fct2019:StationCodeLIS
## Year_fct2016:StationCodeSTTD
## Year_fct2017:StationCodeSTTD
## Year_fct2018:StationCodeSTTD
## Year_fct2019:StationCodeSTTD .
## Year_fct2016:StationCodeLIB
## Year_fct2017:StationCodeLIB
## Year_fct2018:StationCodeLIB **
## Year_fct2019:StationCodeLIB
## Year_fct2016:StationCodeRVB
## Year_fct2017:StationCodeRVB
## Year_fct2018:StationCodeRVB
## Year_fct2019:StationCodeRVB
## Year_fct2016:WaterTemp **
## Year_fct2017:WaterTemp
## Year_fct2018:WaterTemp
## Year_fct2019:WaterTemp
## StationCodeRD22:WaterTemp **
## StationCodeI80:WaterTemp
## StationCodeLIS:WaterTemp *
## StationCodeSTTD:WaterTemp
## StationCodeLIB:WaterTemp *
## StationCodeRVB:WaterTemp
## StationCodeRD22:FlowActionPeriodDuring
## StationCodeI80:FlowActionPeriodDuring
## StationCodeLIS:FlowActionPeriodDuring
## StationCodeSTTD:FlowActionPeriodDuring
## StationCodeLIB:FlowActionPeriodDuring
## StationCodeRVB:FlowActionPeriodDuring
## StationCodeRD22:FlowActionPeriodAfter
## StationCodeI80:FlowActionPeriodAfter
## StationCodeLIS:FlowActionPeriodAfter
## StationCodeSTTD:FlowActionPeriodAfter **
## StationCodeLIB:FlowActionPeriodAfter .
## StationCodeRVB:FlowActionPeriodAfter
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.656 3.656 7.364 1.34e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.762
## Scale est. = 0.32529 n = 3478appraise(m_gamm_cat_wt_ar2b_gam)
k.check(m_gamm_cat_wt_ar2b_gam)## k' edf k-index p-value
## s(DOY) 4 3.655508 0.7254327 0draw(m_gamm_cat_wt_ar2b_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_cat_wt_ar2b_gam, pages = 1, all.terms = TRUE)
acf(resid_cat_wt_ar2b)
Box.test(resid_cat_wt_ar2b, lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: resid_cat_wt_ar2b
## X-squared = 31.046, df = 20, p-value = 0.05459anova(m_gamm_cat_wt_ar2b_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + StationCode + WaterTemp)^2 + StationCode *
## FlowActionPeriod + s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 2.316 0.055024
## StationCode 6 4.867 5.84e-05
## WaterTemp 1 0.526 0.468278
## FlowActionPeriod 2 0.415 0.660676
## Year_fct:StationCode 24 2.036 0.002086
## Year_fct:WaterTemp 4 5.222 0.000342
## StationCode:WaterTemp 6 6.644 5.32e-07
## StationCode:FlowActionPeriod 12 5.012 2.55e-08
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.656 3.656 7.364 1.34e-05This model still looks pretty good. We’ll use this one for the model selection process.
GAM Model 3: Flow as continuous predictor
Plots
df_chla_wt_q %>%
ggplot(aes(x = Flow, y = Chla)) +
geom_point() +
geom_smooth(formula = "y ~ x") +
facet_wrap(vars(StationCode), scales = "free") +
theme_bw()
Let’s break these scatterplots apart by year to see how these patterns vary annually.
df_chla_wt_q %>%
ggplot(aes(x = Flow, y = Chla)) +
geom_point() +
geom_smooth(formula = "y ~ x") +
facet_wrap(
vars(StationCode, Year),
ncol = 5,
scales = "free",
labeller = labeller(.multi_line = FALSE)
) +
theme_bw()
Initial Model
We’ll try running a GAM including all two-way interactions between Year, Daily Average Flow, and Station, and a smooth term for day of year to account for seasonality (restricting the k-value to 5 to reduce overfitting). First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_q <- gam(
Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)Lets look at the model summary and diagnostics:
summary(m_gam_q)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.322e+00 6.781e-02 137.468 < 2e-16 ***
## Year_fct2016 -6.414e-02 7.625e-02 -0.841 0.400279
## Year_fct2017 1.818e-01 7.628e-02 2.383 0.017221 *
## Year_fct2018 1.125e-01 7.096e-02 1.585 0.113065
## Year_fct2019 1.372e-02 7.116e-02 0.193 0.847183
## Flow -1.275e-03 1.059e-04 -12.038 < 2e-16 ***
## StationCodeRD22 3.344e-01 7.804e-02 4.285 1.88e-05 ***
## StationCodeI80 5.559e-01 7.984e-02 6.963 3.98e-12 ***
## StationCodeLIS -6.816e-01 7.585e-02 -8.987 < 2e-16 ***
## StationCodeSTTD -6.447e-01 7.727e-02 -8.344 < 2e-16 ***
## StationCodeLIB -1.640e+00 1.260e-01 -13.019 < 2e-16 ***
## StationCodeRVB -1.908e+00 1.136e-01 -16.801 < 2e-16 ***
## Year_fct2016:Flow -8.355e-05 3.159e-05 -2.645 0.008212 **
## Year_fct2017:Flow -7.085e-05 2.628e-05 -2.696 0.007046 **
## Year_fct2018:Flow -2.494e-05 2.636e-05 -0.946 0.344219
## Year_fct2019:Flow -5.410e-05 2.800e-05 -1.932 0.053386 .
## Year_fct2016:StationCodeRD22 -4.235e-01 9.371e-02 -4.519 6.41e-06 ***
## Year_fct2017:StationCodeRD22 -5.668e-01 9.263e-02 -6.119 1.05e-09 ***
## Year_fct2018:StationCodeRD22 -4.231e-01 8.681e-02 -4.874 1.14e-06 ***
## Year_fct2019:StationCodeRD22 -1.754e-01 8.709e-02 -2.014 0.044067 *
## Year_fct2016:StationCodeI80 2.386e-01 9.525e-02 2.505 0.012300 *
## Year_fct2017:StationCodeI80 -2.309e-01 9.413e-02 -2.453 0.014198 *
## Year_fct2018:StationCodeI80 -4.800e-01 8.848e-02 -5.425 6.21e-08 ***
## Year_fct2019:StationCodeI80 5.216e-01 8.876e-02 5.876 4.60e-09 ***
## Year_fct2016:StationCodeLIS 5.173e-01 9.177e-02 5.637 1.87e-08 ***
## Year_fct2017:StationCodeLIS 3.056e-01 9.172e-02 3.332 0.000870 ***
## Year_fct2018:StationCodeLIS -1.288e-01 8.627e-02 -1.493 0.135495
## Year_fct2019:StationCodeLIS 5.536e-01 8.765e-02 6.316 3.02e-10 ***
## Year_fct2016:StationCodeSTTD 3.119e-01 9.445e-02 3.302 0.000971 ***
## Year_fct2017:StationCodeSTTD -1.372e-01 9.943e-02 -1.380 0.167655
## Year_fct2018:StationCodeSTTD -7.586e-01 9.080e-02 -8.354 < 2e-16 ***
## Year_fct2019:StationCodeSTTD -1.145e+00 8.934e-02 -12.815 < 2e-16 ***
## Year_fct2016:StationCodeLIB 8.329e-01 1.058e-01 7.871 4.67e-15 ***
## Year_fct2017:StationCodeLIB -1.933e-01 1.063e-01 -1.818 0.069088 .
## Year_fct2018:StationCodeLIB -1.979e+00 1.223e-01 -16.186 < 2e-16 ***
## Year_fct2019:StationCodeLIB -7.144e-01 1.241e-01 -5.757 9.31e-09 ***
## Year_fct2016:StationCodeRVB 6.335e-01 2.347e-01 2.699 0.006986 **
## Year_fct2017:StationCodeRVB -2.198e-02 1.987e-01 -0.111 0.911935
## Year_fct2018:StationCodeRVB -4.140e-01 1.626e-01 -2.546 0.010955 *
## Year_fct2019:StationCodeRVB -4.845e-01 1.910e-01 -2.537 0.011238 *
## Flow:StationCodeRD22 -1.057e-04 1.375e-04 -0.768 0.442446
## Flow:StationCodeI80 -6.752e-04 1.387e-04 -4.867 1.18e-06 ***
## Flow:StationCodeLIS 1.541e-03 1.391e-04 11.079 < 2e-16 ***
## Flow:StationCodeSTTD 3.211e-03 1.406e-04 22.832 < 2e-16 ***
## Flow:StationCodeLIB 1.275e-03 1.220e-04 10.452 < 2e-16 ***
## Flow:StationCodeRVB 1.327e-03 1.039e-04 12.773 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.757 3.966 81.39 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.851 Deviance explained = 85.3%
## -REML = 1770.4 Scale est. = 0.14725 n = 3478appraise(m_gam_q)
k.check(m_gam_q)## k' edf k-index p-value
## s(DOY) 4 3.756521 0.9228792 0draw(m_gam_q, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_q, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_q))
Box.test(residuals(m_gam_q), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_q)
## X-squared = 10096, df = 20, p-value < 2.2e-16Model with AR() correlation structure
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_q <- as.formula("Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)")
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q <- gamm(
f_gam_q,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_q_ar1 <- gamm(
f_gam_q,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_ar2 <- gamm(
f_gam_q,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_ar3 <- gamm(
f_gam_q,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_q$lme,
m_gamm_q_ar1$lme,
m_gamm_q_ar2$lme,
m_gamm_q_ar3$lme
)## Model df AIC BIC logLik Test L.Ratio p-value
## m_gamm_q$lme 1 49 3638.827 3939.716 -1770.413
## m_gamm_q_ar1$lme 2 50 -2092.217 -1785.187 1096.109 1 vs 2 5733.044 <.0001
## m_gamm_q_ar2$lme 3 51 -2111.366 -1798.195 1106.683 2 vs 3 21.149 <.0001
## m_gamm_q_ar3$lme 4 52 -2110.010 -1790.698 1107.005 3 vs 4 0.644 0.4224It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
AR(2) Model
# Pull out the residuals and the GAM model
resid_q_ar2 <- residuals(m_gamm_q_ar2$lme, type = "normalized")
m_gamm_q_ar2_gam <- m_gamm_q_ar2$gam
summary(m_gamm_q_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.178e+00 2.921e-01 31.417 < 2e-16 ***
## Year_fct2016 -3.270e-02 3.761e-01 -0.087 0.930719
## Year_fct2017 3.329e-01 3.587e-01 0.928 0.353319
## Year_fct2018 1.306e-01 3.495e-01 0.374 0.708602
## Year_fct2019 8.786e-02 3.505e-01 0.251 0.802084
## Flow -1.821e-04 1.772e-04 -1.028 0.304126
## StationCodeRD22 3.170e-01 3.570e-01 0.888 0.374697
## StationCodeI80 7.103e-01 3.671e-01 1.935 0.053071 .
## StationCodeLIS -5.832e-01 3.496e-01 -1.668 0.095356 .
## StationCodeSTTD -4.766e-01 3.576e-01 -1.333 0.182705
## StationCodeLIB -1.228e+00 3.569e-01 -3.440 0.000588 ***
## StationCodeRVB -1.517e+00 3.519e-01 -4.312 1.66e-05 ***
## Year_fct2016:Flow -5.265e-06 1.924e-05 -0.274 0.784344
## Year_fct2017:Flow -4.817e-06 1.964e-05 -0.245 0.806252
## Year_fct2018:Flow -3.172e-06 2.102e-05 -0.151 0.880067
## Year_fct2019:Flow -2.436e-05 1.961e-05 -1.242 0.214151
## Year_fct2016:StationCodeRD22 -4.412e-01 4.823e-01 -0.915 0.360450
## Year_fct2017:StationCodeRD22 -5.891e-01 4.591e-01 -1.283 0.199529
## Year_fct2018:StationCodeRD22 -4.262e-01 4.486e-01 -0.950 0.342068
## Year_fct2019:StationCodeRD22 -1.675e-01 4.503e-01 -0.372 0.709983
## Year_fct2016:StationCodeI80 1.284e-01 4.896e-01 0.262 0.793087
## Year_fct2017:StationCodeI80 -3.791e-01 4.669e-01 -0.812 0.416906
## Year_fct2018:StationCodeI80 -6.219e-01 4.564e-01 -1.363 0.173092
## Year_fct2019:StationCodeI80 3.949e-01 4.582e-01 0.862 0.388749
## Year_fct2016:StationCodeLIS 6.885e-01 4.700e-01 1.465 0.143014
## Year_fct2017:StationCodeLIS 1.647e-01 4.566e-01 0.361 0.718331
## Year_fct2018:StationCodeLIS -1.577e-01 4.445e-01 -0.355 0.722725
## Year_fct2019:StationCodeLIS 4.135e-01 4.496e-01 0.920 0.357786
## Year_fct2016:StationCodeSTTD 3.444e-01 4.849e-01 0.710 0.477618
## Year_fct2017:StationCodeSTTD -5.748e-01 4.914e-01 -1.170 0.242151
## Year_fct2018:StationCodeSTTD -8.724e-01 4.660e-01 -1.872 0.061252 .
## Year_fct2019:StationCodeSTTD -1.196e+00 4.584e-01 -2.609 0.009123 **
## Year_fct2016:StationCodeLIB 9.423e-01 4.738e-01 1.989 0.046817 *
## Year_fct2017:StationCodeLIB -3.952e-01 4.620e-01 -0.855 0.392434
## Year_fct2018:StationCodeLIB -2.200e+00 5.051e-01 -4.356 1.36e-05 ***
## Year_fct2019:StationCodeLIB -9.373e-01 5.127e-01 -1.828 0.067617 .
## Year_fct2016:StationCodeRVB 3.293e-01 5.016e-01 0.657 0.511538
## Year_fct2017:StationCodeRVB -5.462e-01 4.829e-01 -1.131 0.258170
## Year_fct2018:StationCodeRVB -4.555e-01 4.624e-01 -0.985 0.324691
## Year_fct2019:StationCodeRVB -6.753e-01 4.697e-01 -1.438 0.150619
## Flow:StationCodeRD22 -4.185e-04 2.542e-04 -1.646 0.099770 .
## Flow:StationCodeI80 -1.249e-03 2.832e-04 -4.412 1.06e-05 ***
## Flow:StationCodeLIS -2.168e-05 2.860e-04 -0.076 0.939582
## Flow:StationCodeSTTD 1.431e-03 2.903e-04 4.930 8.61e-07 ***
## Flow:StationCodeLIB 3.575e-04 1.821e-04 1.964 0.049625 *
## Flow:StationCodeRVB 1.822e-04 1.771e-04 1.029 0.303658
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.675 3.675 11.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.819
## Scale est. = 0.22342 n = 3478appraise(m_gamm_q_ar2_gam)
k.check(m_gamm_q_ar2_gam)## k' edf k-index p-value
## s(DOY) 4 3.674585 0.803638 0draw(m_gamm_q_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_q_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_q_ar2)
Box.test(resid_q_ar2, lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: resid_q_ar2
## X-squared = 29.682, df = 20, p-value = 0.07517The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 0.407 0.803
## Flow 1 1.056 0.304
## StationCode 6 14.929 < 2e-16
## Year_fct:Flow 4 0.408 0.803
## Year_fct:StationCode 24 3.933 3.9e-10
## Flow:StationCode 6 16.764 < 2e-16
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.675 3.675 11.3 <2e-16The Year:StationCode and Flow:StationCode interactions were significant among the interaction terms of the model. Let’s remove the Year:Flow interaction term and check the model summary and diagnostics.
m_gamm_q_ar2b <- gamm(
Chla_log ~ Year_fct * StationCode + Flow * StationCode + s(DOY, k = 5),
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
# Pull out the residuals and the GAM model
resid_q_ar2b <- residuals(m_gamm_q_ar2b$lme, type = "normalized")
m_gamm_q_ar2b_gam <- m_gamm_q_ar2b$gam
summary(m_gamm_q_ar2b_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + s(DOY,
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.180e+00 2.920e-01 31.445 < 2e-16 ***
## Year_fct2016 -3.484e-02 3.758e-01 -0.093 0.926138
## Year_fct2017 3.299e-01 3.584e-01 0.920 0.357500
## Year_fct2018 1.278e-01 3.493e-01 0.366 0.714492
## Year_fct2019 8.237e-02 3.503e-01 0.235 0.814093
## StationCodeRD22 3.144e-01 3.568e-01 0.881 0.378251
## StationCodeI80 7.079e-01 3.669e-01 1.930 0.053736 .
## StationCodeLIS -5.870e-01 3.494e-01 -1.680 0.093030 .
## StationCodeSTTD -4.804e-01 3.574e-01 -1.344 0.178970
## StationCodeLIB -1.239e+00 3.563e-01 -3.477 0.000514 ***
## StationCodeRVB -1.490e+00 3.492e-01 -4.266 2.04e-05 ***
## Flow -1.903e-04 1.769e-04 -1.076 0.282092
## Year_fct2016:StationCodeRD22 -4.388e-01 4.821e-01 -0.910 0.362757
## Year_fct2017:StationCodeRD22 -5.866e-01 4.589e-01 -1.278 0.201177
## Year_fct2018:StationCodeRD22 -4.239e-01 4.483e-01 -0.946 0.344444
## Year_fct2019:StationCodeRD22 -1.649e-01 4.501e-01 -0.366 0.714156
## Year_fct2016:StationCodeI80 1.310e-01 4.893e-01 0.268 0.788936
## Year_fct2017:StationCodeI80 -3.767e-01 4.666e-01 -0.807 0.419558
## Year_fct2018:StationCodeI80 -6.194e-01 4.562e-01 -1.358 0.174616
## Year_fct2019:StationCodeI80 3.978e-01 4.579e-01 0.869 0.385022
## Year_fct2016:StationCodeLIS 6.924e-01 4.697e-01 1.474 0.140536
## Year_fct2017:StationCodeLIS 1.681e-01 4.563e-01 0.368 0.712634
## Year_fct2018:StationCodeLIS -1.545e-01 4.442e-01 -0.348 0.728056
## Year_fct2019:StationCodeLIS 4.184e-01 4.493e-01 0.931 0.351866
## Year_fct2016:StationCodeSTTD 3.480e-01 4.846e-01 0.718 0.472790
## Year_fct2017:StationCodeSTTD -5.708e-01 4.911e-01 -1.162 0.245159
## Year_fct2018:StationCodeSTTD -8.680e-01 4.657e-01 -1.864 0.062400 .
## Year_fct2019:StationCodeSTTD -1.191e+00 4.581e-01 -2.600 0.009364 **
## Year_fct2016:StationCodeLIB 9.536e-01 4.727e-01 2.017 0.043741 *
## Year_fct2017:StationCodeLIB -3.846e-01 4.611e-01 -0.834 0.404321
## Year_fct2018:StationCodeLIB -2.191e+00 5.042e-01 -4.345 1.43e-05 ***
## Year_fct2019:StationCodeLIB -9.096e-01 5.118e-01 -1.777 0.075641 .
## Year_fct2016:StationCodeRVB 3.114e-01 4.819e-01 0.646 0.518156
## Year_fct2017:StationCodeRVB -5.575e-01 4.532e-01 -1.230 0.218786
## Year_fct2018:StationCodeRVB -4.597e-01 4.419e-01 -1.040 0.298248
## Year_fct2019:StationCodeRVB -8.461e-01 4.498e-01 -1.881 0.060059 .
## StationCodeRD22:Flow -4.181e-04 2.541e-04 -1.645 0.100007
## StationCodeI80:Flow -1.252e-03 2.831e-04 -4.424 1.00e-05 ***
## StationCodeLIS:Flow -2.371e-05 2.859e-04 -0.083 0.933893
## StationCodeSTTD:Flow 1.429e-03 2.902e-04 4.924 8.88e-07 ***
## StationCodeLIB:Flow 3.608e-04 1.819e-04 1.983 0.047421 *
## StationCodeRVB:Flow 1.840e-04 1.770e-04 1.039 0.298683
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.68 3.68 11.47 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.819
## Scale est. = 0.22326 n = 3478appraise(m_gamm_q_ar2b_gam)
k.check(m_gamm_q_ar2b_gam)## k' edf k-index p-value
## s(DOY) 4 3.679956 0.8026532 0draw(m_gamm_q_ar2b_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_q_ar2b_gam, pages = 1, all.terms = TRUE)
acf(resid_q_ar2b)
Box.test(resid_q_ar2b, lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: resid_q_ar2b
## X-squared = 29.646, df = 20, p-value = 0.0758anova(m_gamm_q_ar2b_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + s(DOY,
## k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 0.406 0.805
## StationCode 6 14.873 < 2e-16
## Flow 1 1.157 0.282
## Year_fct:StationCode 24 4.102 8.4e-11
## StationCode:Flow 6 16.865 < 2e-16
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.68 3.68 11.47 <2e-16This model still looks pretty good. We’ll use this one for the model selection process.
GAM Model 4: Flow and water temperature as continuous predictors
Plots
Before we add both flow and water temperature as continuous predictors in the model, let’s check if they are correlated with each other.
df_chla_wt_q %>%
ggplot(aes(x = Flow, y = WaterTemp)) +
geom_point() +
geom_smooth(formula = "y ~ x") +
facet_wrap(vars(StationCode), scales = "free") +
theme_bw()
Flow and water temperature appear to have no correlation.
Initial Model
We’ll try running a GAM using the model structure from the flow model and adding two-way interactions between Year:WaterTemp and Station:WaterTemp. First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_q_wt <- gam(
Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)Lets look at the model summary and diagnostics:
summary(m_gam_q_wt)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct *
## WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.2662032 0.2178374 47.128 < 2e-16 ***
## Year_fct2016 -0.6484169 0.3059264 -2.120 0.034118 *
## Year_fct2017 -0.0609459 0.1748727 -0.349 0.727474
## Year_fct2018 -0.1613266 0.1680842 -0.960 0.337226
## Year_fct2019 0.5338094 0.1622929 3.289 0.001015 **
## StationCodeRD22 -1.0478545 0.1618085 -6.476 1.08e-10 ***
## StationCodeI80 -0.3198595 0.1552094 -2.061 0.039395 *
## StationCodeLIS -2.5319716 0.1616506 -15.663 < 2e-16 ***
## StationCodeSTTD -1.7388040 0.1881330 -9.242 < 2e-16 ***
## StationCodeLIB -3.9184097 0.2367461 -16.551 < 2e-16 ***
## StationCodeRVB -3.4749646 0.1941132 -17.902 < 2e-16 ***
## Flow -0.0011256 0.0001014 -11.096 < 2e-16 ***
## WaterTemp -0.0433571 0.0092626 -4.681 2.97e-06 ***
## Year_fct2016:StationCodeRD22 -0.5347596 0.0910571 -5.873 4.69e-09 ***
## Year_fct2017:StationCodeRD22 -0.6200951 0.0889603 -6.970 3.77e-12 ***
## Year_fct2018:StationCodeRD22 -0.3517822 0.0834277 -4.217 2.54e-05 ***
## Year_fct2019:StationCodeRD22 -0.1482813 0.0833741 -1.779 0.075410 .
## Year_fct2016:StationCodeI80 0.2012539 0.0934108 2.155 0.031270 *
## Year_fct2017:StationCodeI80 -0.2400564 0.0912588 -2.631 0.008564 **
## Year_fct2018:StationCodeI80 -0.4104957 0.0855122 -4.800 1.65e-06 ***
## Year_fct2019:StationCodeI80 0.5658673 0.0856073 6.610 4.44e-11 ***
## Year_fct2016:StationCodeLIS 0.3694915 0.0892519 4.140 3.56e-05 ***
## Year_fct2017:StationCodeLIS 0.2490534 0.0878966 2.833 0.004631 **
## Year_fct2018:StationCodeLIS -0.0695068 0.0825594 -0.842 0.399902
## Year_fct2019:StationCodeLIS 0.5781116 0.0838872 6.892 6.54e-12 ***
## Year_fct2016:StationCodeSTTD 0.2714631 0.0931049 2.916 0.003572 **
## Year_fct2017:StationCodeSTTD -0.1673873 0.0975732 -1.716 0.086343 .
## Year_fct2018:StationCodeSTTD -0.6929241 0.0871555 -7.950 2.50e-15 ***
## Year_fct2019:StationCodeSTTD -1.1345661 0.0858236 -13.220 < 2e-16 ***
## Year_fct2016:StationCodeLIB 0.8903568 0.0949408 9.378 < 2e-16 ***
## Year_fct2017:StationCodeLIB -0.2315049 0.0969164 -2.389 0.016962 *
## Year_fct2018:StationCodeLIB -2.0080384 0.1150845 -17.448 < 2e-16 ***
## Year_fct2019:StationCodeLIB -0.8582224 0.1196082 -7.175 8.81e-13 ***
## Year_fct2016:StationCodeRVB 0.1597391 0.1027829 1.554 0.120243
## Year_fct2017:StationCodeRVB -0.3753724 0.1038743 -3.614 0.000306 ***
## Year_fct2018:StationCodeRVB -0.2920268 0.0889282 -3.284 0.001034 **
## Year_fct2019:StationCodeRVB -0.6540161 0.0928199 -7.046 2.22e-12 ***
## StationCodeRD22:Flow -0.0003430 0.0001340 -2.560 0.010519 *
## StationCodeI80:Flow -0.0008049 0.0001352 -5.951 2.93e-09 ***
## StationCodeLIS:Flow 0.0012604 0.0001350 9.336 < 2e-16 ***
## StationCodeSTTD:Flow 0.0030649 0.0001386 22.121 < 2e-16 ***
## StationCodeLIB:Flow 0.0013482 0.0001236 10.908 < 2e-16 ***
## StationCodeRVB:Flow 0.0011204 0.0001014 11.046 < 2e-16 ***
## Year_fct2016:WaterTemp 0.0278927 0.0123785 2.253 0.024302 *
## Year_fct2017:WaterTemp 0.0125116 0.0069925 1.789 0.073658 .
## Year_fct2018:WaterTemp 0.0103364 0.0068995 1.498 0.134191
## Year_fct2019:WaterTemp -0.0253859 0.0065403 -3.881 0.000106 ***
## StationCodeRD22:WaterTemp 0.0637778 0.0067035 9.514 < 2e-16 ***
## StationCodeI80:WaterTemp 0.0389816 0.0064179 6.074 1.39e-09 ***
## StationCodeLIS:WaterTemp 0.0857066 0.0066872 12.816 < 2e-16 ***
## StationCodeSTTD:WaterTemp 0.0496447 0.0081526 6.089 1.26e-09 ***
## StationCodeLIB:WaterTemp 0.1230615 0.0120862 10.182 < 2e-16 ***
## StationCodeRVB:WaterTemp 0.0833896 0.0086653 9.623 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.781 3.974 53.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.864 Deviance explained = 86.6%
## -REML = 1615.8 Scale est. = 0.13435 n = 3478appraise(m_gam_q_wt)
k.check(m_gam_q_wt)## k' edf k-index p-value
## s(DOY) 4 3.78136 0.9264289 0draw(m_gam_q_wt, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_q_wt, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_q_wt))
Box.test(residuals(m_gam_q_wt), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_q_wt)
## X-squared = 8803.9, df = 20, p-value < 2.2e-16Model with AR() correlation structure
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_q_wt <- as.formula(
"Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)"
)
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q_wt <- gamm(
f_gam_q_wt,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_q_wt_ar1 <- gamm(
f_gam_q_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_wt_ar2 <- gamm(
f_gam_q_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_wt_ar3 <- gamm(
f_gam_q_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_q_wt$lme,
m_gamm_q_wt_ar1$lme,
m_gamm_q_wt_ar2$lme,
m_gamm_q_wt_ar3$lme
)## Model df AIC BIC logLik Test L.Ratio
## m_gamm_q_wt$lme 1 56 3343.533 3687.292 -1615.766
## m_gamm_q_wt_ar1$lme 2 57 -2152.733 -1802.835 1133.367 1 vs 2 5498.266
## m_gamm_q_wt_ar2$lme 3 58 -2174.391 -1818.354 1145.195 2 vs 3 23.657
## m_gamm_q_wt_ar3$lme 4 59 -2173.370 -1811.195 1145.685 3 vs 4 0.980
## p-value
## m_gamm_q_wt$lme
## m_gamm_q_wt_ar1$lme <.0001
## m_gamm_q_wt_ar2$lme <.0001
## m_gamm_q_wt_ar3$lme 0.3222It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
AR(2) Model
# Pull out the residuals and the GAM model
resid_q_wt_ar2 <- residuals(m_gamm_q_wt_ar2$lme, type = "normalized")
m_gamm_q_wt_ar2_gam <- m_gamm_q_wt_ar2$gam
summary(m_gamm_q_wt_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct *
## WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.894e+00 5.959e-01 14.925 < 2e-16 ***
## Year_fct2016 1.423e+00 6.306e-01 2.256 0.024113 *
## Year_fct2017 4.921e-02 5.714e-01 0.086 0.931374
## Year_fct2018 -3.647e-01 5.778e-01 -0.631 0.527977
## Year_fct2019 -2.601e-02 5.490e-01 -0.047 0.962218
## StationCodeRD22 -1.128e+00 6.141e-01 -1.838 0.066189 .
## StationCodeI80 8.206e-01 5.957e-01 1.378 0.168447
## StationCodeLIS -1.605e+00 6.037e-01 -2.659 0.007879 **
## StationCodeSTTD -1.252e+00 6.252e-01 -2.003 0.045228 *
## StationCodeLIB -9.809e-01 6.957e-01 -1.410 0.158668
## StationCodeRVB -2.142e+00 7.531e-01 -2.845 0.004470 **
## Flow -1.699e-04 1.767e-04 -0.961 0.336414
## WaterTemp 1.399e-02 2.161e-02 0.648 0.517248
## Year_fct2016:StationCodeRD22 -6.003e-01 5.380e-01 -1.116 0.264600
## Year_fct2017:StationCodeRD22 -6.123e-01 5.116e-01 -1.197 0.231484
## Year_fct2018:StationCodeRD22 -3.207e-01 5.017e-01 -0.639 0.522705
## Year_fct2019:StationCodeRD22 -1.397e-01 5.030e-01 -0.278 0.781210
## Year_fct2016:StationCodeI80 2.047e-01 5.478e-01 0.374 0.708640
## Year_fct2017:StationCodeI80 -2.886e-01 5.216e-01 -0.553 0.580145
## Year_fct2018:StationCodeI80 -5.347e-01 5.115e-01 -1.045 0.295920
## Year_fct2019:StationCodeI80 4.585e-01 5.129e-01 0.894 0.371361
## Year_fct2016:StationCodeLIS 5.821e-01 5.244e-01 1.110 0.267053
## Year_fct2017:StationCodeLIS 1.927e-01 5.089e-01 0.379 0.704890
## Year_fct2018:StationCodeLIS -1.008e-01 4.971e-01 -0.203 0.839336
## Year_fct2019:StationCodeLIS 4.263e-01 5.023e-01 0.849 0.396085
## Year_fct2016:StationCodeSTTD 2.585e-01 5.418e-01 0.477 0.633378
## Year_fct2017:StationCodeSTTD -5.788e-01 5.483e-01 -1.056 0.291233
## Year_fct2018:StationCodeSTTD -8.371e-01 5.208e-01 -1.607 0.108099
## Year_fct2019:StationCodeSTTD -1.172e+00 5.124e-01 -2.288 0.022224 *
## Year_fct2016:StationCodeLIB 8.689e-01 5.287e-01 1.644 0.100355
## Year_fct2017:StationCodeLIB -2.787e-01 5.147e-01 -0.541 0.588263
## Year_fct2018:StationCodeLIB -2.096e+00 5.618e-01 -3.731 0.000194 ***
## Year_fct2019:StationCodeLIB -8.120e-01 5.685e-01 -1.428 0.153286
## Year_fct2016:StationCodeRVB 1.962e-01 5.384e-01 0.364 0.715519
## Year_fct2017:StationCodeRVB -4.613e-01 5.059e-01 -0.912 0.361903
## Year_fct2018:StationCodeRVB -3.546e-01 4.956e-01 -0.715 0.474415
## Year_fct2019:StationCodeRVB -7.574e-01 5.041e-01 -1.503 0.133011
## StationCodeRD22:Flow -4.412e-04 2.550e-04 -1.730 0.083665 .
## StationCodeI80:Flow -1.205e-03 2.852e-04 -4.226 2.44e-05 ***
## StationCodeLIS:Flow -8.722e-06 2.877e-04 -0.030 0.975817
## StationCodeSTTD:Flow 1.426e-03 2.923e-04 4.878 1.12e-06 ***
## StationCodeLIB:Flow 3.233e-04 1.821e-04 1.775 0.075912 .
## StationCodeRVB:Flow 1.655e-04 1.768e-04 0.936 0.349418
## Year_fct2016:WaterTemp -5.896e-02 2.002e-02 -2.945 0.003250 **
## Year_fct2017:WaterTemp 1.079e-02 1.791e-02 0.603 0.546631
## Year_fct2018:WaterTemp 2.128e-02 1.913e-02 1.113 0.265978
## Year_fct2019:WaterTemp 3.786e-03 1.708e-02 0.222 0.824591
## StationCodeRD22:WaterTemp 6.651e-02 2.099e-02 3.169 0.001541 **
## StationCodeI80:WaterTemp -8.905e-03 1.954e-02 -0.456 0.648553
## StationCodeLIS:WaterTemp 4.642e-02 2.053e-02 2.261 0.023804 *
## StationCodeSTTD:WaterTemp 3.510e-02 2.186e-02 1.606 0.108445
## StationCodeLIB:WaterTemp -1.635e-02 2.705e-02 -0.605 0.545485
## StationCodeRVB:WaterTemp 2.870e-02 2.960e-02 0.970 0.332241
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.664 3.664 8.628 2.2e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.807
## Scale est. = 0.2515 n = 3478appraise(m_gamm_q_wt_ar2_gam)
k.check(m_gamm_q_wt_ar2_gam)## k' edf k-index p-value
## s(DOY) 4 3.663875 0.7855439 0draw(m_gamm_q_wt_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_q_wt_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_q_wt_ar2)
Box.test(resid_q_wt_ar2, lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: resid_q_wt_ar2
## X-squared = 30.519, df = 20, p-value = 0.06187The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_wt_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct *
## WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 2.531 0.038608
## StationCode 6 6.652 5.20e-07
## Flow 1 0.924 0.336414
## WaterTemp 1 0.419 0.517248
## Year_fct:StationCode 24 3.165 3.27e-07
## StationCode:Flow 6 15.268 < 2e-16
## Year_fct:WaterTemp 4 4.659 0.000945
## StationCode:WaterTemp 6 4.669 9.75e-05
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.664 3.664 8.628 2.2e-06All interaction terms are significant in this model, so we’ll use this one for the model selection process.
GAM Model 5: Flow with flow action period
Plots
Let’s look at the relationship between chlorophyll and flow for each station-flow action period combination.
df_chla_wt_q %>%
ggplot(aes(x = Flow, y = Chla)) +
geom_point() +
geom_smooth(formula = "y ~ x") +
facet_wrap(
vars(StationCode, FlowActionPeriod),
ncol = 3,
scales = "free",
labeller = labeller(.multi_line = FALSE)
) +
theme_bw()
Initial Model
It’s possible that the relationship between chlorophyll and Flow differs between Flow Action Periods at each station. We’ll try to account for this by running a GAM including a two-way interaction between Year and Station, a three-way interaction between Flow, Station, and Flow Action Period, and a smooth term for day of year to account for seasonality (restricting the k-value to 5 to reduce overfitting). We won’t include the interaction between Year and Flow Action Period since it was barely significant in the categorical model and we’re primarily including the Flow Action Period to allow for different slopes for flow at each station. First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_q_fa <- gam(
Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)Lets look at the model summary and diagnostics:
summary(m_gam_q_fa)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value
## (Intercept) 9.0200439 0.0839467 107.450
## Year_fct2016 -0.0752673 0.0817958 -0.920
## Year_fct2017 0.1772134 0.0707049 2.506
## Year_fct2018 -0.0051101 0.0673845 -0.076
## Year_fct2019 -0.0504808 0.0701749 -0.719
## StationCodeRD22 0.6048920 0.1021743 5.920
## StationCodeI80 0.7075466 0.1053335 6.717
## StationCodeLIS -0.2534964 0.1008109 -2.515
## StationCodeSTTD -0.4888704 0.1077499 -4.537
## StationCodeLIB -0.5233575 0.1654118 -3.164
## StationCodeRVB -1.3296625 0.1140411 -11.660
## Flow 0.0003845 0.0006336 0.607
## FlowActionPeriodDuring 0.2206777 0.0876376 2.518
## FlowActionPeriodAfter 0.8225799 0.0753879 10.911
## Year_fct2016:StationCodeRD22 -0.3499073 0.0955876 -3.661
## Year_fct2017:StationCodeRD22 -0.5693333 0.0856991 -6.643
## Year_fct2018:StationCodeRD22 -0.3195979 0.0814168 -3.925
## Year_fct2019:StationCodeRD22 -0.0887040 0.0850886 -1.042
## Year_fct2016:StationCodeI80 0.3049614 0.0984091 3.099
## Year_fct2017:StationCodeI80 -0.2386526 0.0866819 -2.753
## Year_fct2018:StationCodeI80 -0.3895701 0.0840942 -4.633
## Year_fct2019:StationCodeI80 0.6211290 0.0862887 7.198
## Year_fct2016:StationCodeLIS 0.3767691 0.0932856 4.039
## Year_fct2017:StationCodeLIS 0.2509619 0.0859328 2.920
## Year_fct2018:StationCodeLIS 0.0062195 0.0810131 0.077
## Year_fct2019:StationCodeLIS 0.6367041 0.0862458 7.382
## Year_fct2016:StationCodeSTTD 0.2375629 0.0960514 2.473
## Year_fct2017:StationCodeSTTD -0.0453964 0.0976543 -0.465
## Year_fct2018:StationCodeSTTD -0.5512862 0.0857425 -6.430
## Year_fct2019:StationCodeSTTD -1.0150911 0.0885178 -11.468
## Year_fct2016:StationCodeLIB 0.8861875 0.0923868 9.592
## Year_fct2017:StationCodeLIB -0.2709091 0.0930803 -2.910
## Year_fct2018:StationCodeLIB -2.0402939 0.1103463 -18.490
## Year_fct2019:StationCodeLIB -1.0231592 0.1184799 -8.636
## Year_fct2016:StationCodeRVB 0.1576659 0.1008718 1.563
## Year_fct2017:StationCodeRVB -0.4120981 0.0981257 -4.200
## Year_fct2018:StationCodeRVB -0.2930854 0.0862916 -3.396
## Year_fct2019:StationCodeRVB -0.6454502 0.0914059 -7.061
## StationCodeRD22:Flow 0.0003861 0.0014058 0.275
## StationCodeI80:Flow 0.0018616 0.0018957 0.982
## StationCodeLIS:Flow 0.0030944 0.0009665 3.202
## StationCodeSTTD:Flow 0.0023417 0.0011000 2.129
## StationCodeLIB:Flow 0.0001091 0.0006393 0.171
## StationCodeRVB:Flow -0.0003875 0.0006337 -0.612
## Flow:FlowActionPeriodDuring -0.0013504 0.0006478 -2.084
## Flow:FlowActionPeriodAfter -0.0083453 0.0011337 -7.361
## StationCodeRD22:FlowActionPeriodDuring -0.2621276 0.1125150 -2.330
## StationCodeI80:FlowActionPeriodDuring -0.1046615 0.1170017 -0.895
## StationCodeLIS:FlowActionPeriodDuring 0.3486851 0.1118607 3.117
## StationCodeSTTD:FlowActionPeriodDuring 0.4101368 0.1221588 3.357
## StationCodeLIB:FlowActionPeriodDuring -0.3116617 0.1631523 -1.910
## StationCodeRVB:FlowActionPeriodDuring -0.2557270 0.1403998 -1.821
## StationCodeRD22:FlowActionPeriodAfter -0.6738249 0.0876828 -7.685
## StationCodeI80:FlowActionPeriodAfter -0.4696770 0.0926079 -5.072
## StationCodeLIS:FlowActionPeriodAfter -0.6118818 0.0862036 -7.098
## StationCodeSTTD:FlowActionPeriodAfter -0.4414052 0.0972779 -4.538
## StationCodeLIB:FlowActionPeriodAfter -1.9012296 0.1443312 -13.173
## StationCodeRVB:FlowActionPeriodAfter -0.9357160 0.1208432 -7.743
## StationCodeRD22:Flow:FlowActionPeriodDuring -0.0007004 0.0014165 -0.494
## StationCodeI80:Flow:FlowActionPeriodDuring -0.0028106 0.0019000 -1.479
## StationCodeLIS:Flow:FlowActionPeriodDuring -0.0035976 0.0009869 -3.645
## StationCodeSTTD:Flow:FlowActionPeriodDuring -0.0007149 0.0011208 -0.638
## StationCodeLIB:Flow:FlowActionPeriodDuring 0.0011745 0.0006559 1.791
## StationCodeRVB:Flow:FlowActionPeriodDuring 0.0013536 0.0006479 2.089
## StationCodeRD22:Flow:FlowActionPeriodAfter -0.0017581 0.0023433 -0.750
## StationCodeI80:Flow:FlowActionPeriodAfter -0.0055463 0.0025798 -2.150
## StationCodeLIS:Flow:FlowActionPeriodAfter 0.0125081 0.0015103 8.282
## StationCodeSTTD:Flow:FlowActionPeriodAfter 0.0118661 0.0016282 7.288
## StationCodeLIB:Flow:FlowActionPeriodAfter 0.0076975 0.0011391 6.758
## StationCodeRVB:Flow:FlowActionPeriodAfter 0.0083444 0.0011336 7.361
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year_fct2016 0.357541
## Year_fct2017 0.012244 *
## Year_fct2018 0.939555
## Year_fct2019 0.471970
## StationCodeRD22 3.53e-09 ***
## StationCodeI80 2.16e-11 ***
## StationCodeLIS 0.011963 *
## StationCodeSTTD 5.90e-06 ***
## StationCodeLIB 0.001570 **
## StationCodeRVB < 2e-16 ***
## Flow 0.543936
## FlowActionPeriodDuring 0.011846 *
## FlowActionPeriodAfter < 2e-16 ***
## Year_fct2016:StationCodeRD22 0.000255 ***
## Year_fct2017:StationCodeRD22 3.55e-11 ***
## Year_fct2018:StationCodeRD22 8.83e-05 ***
## Year_fct2019:StationCodeRD22 0.297259
## Year_fct2016:StationCodeI80 0.001958 **
## Year_fct2017:StationCodeI80 0.005933 **
## Year_fct2018:StationCodeI80 3.75e-06 ***
## Year_fct2019:StationCodeI80 7.47e-13 ***
## Year_fct2016:StationCodeLIS 5.49e-05 ***
## Year_fct2017:StationCodeLIS 0.003518 **
## Year_fct2018:StationCodeLIS 0.938810
## Year_fct2019:StationCodeLIS 1.94e-13 ***
## Year_fct2016:StationCodeSTTD 0.013436 *
## Year_fct2017:StationCodeSTTD 0.642055
## Year_fct2018:StationCodeSTTD 1.46e-10 ***
## Year_fct2019:StationCodeSTTD < 2e-16 ***
## Year_fct2016:StationCodeLIB < 2e-16 ***
## Year_fct2017:StationCodeLIB 0.003632 **
## Year_fct2018:StationCodeLIB < 2e-16 ***
## Year_fct2019:StationCodeLIB < 2e-16 ***
## Year_fct2016:StationCodeRVB 0.118138
## Year_fct2017:StationCodeRVB 2.74e-05 ***
## Year_fct2018:StationCodeRVB 0.000690 ***
## Year_fct2019:StationCodeRVB 1.99e-12 ***
## StationCodeRD22:Flow 0.783602
## StationCodeI80:Flow 0.326154
## StationCodeLIS:Flow 0.001378 **
## StationCodeSTTD:Flow 0.033342 *
## StationCodeLIB:Flow 0.864475
## StationCodeRVB:Flow 0.540909
## Flow:FlowActionPeriodDuring 0.037198 *
## Flow:FlowActionPeriodAfter 2.27e-13 ***
## StationCodeRD22:FlowActionPeriodDuring 0.019880 *
## StationCodeI80:FlowActionPeriodDuring 0.371102
## StationCodeLIS:FlowActionPeriodDuring 0.001841 **
## StationCodeSTTD:FlowActionPeriodDuring 0.000795 ***
## StationCodeLIB:FlowActionPeriodDuring 0.056185 .
## StationCodeRVB:FlowActionPeriodDuring 0.068631 .
## StationCodeRD22:FlowActionPeriodAfter 1.99e-14 ***
## StationCodeI80:FlowActionPeriodAfter 4.15e-07 ***
## StationCodeLIS:FlowActionPeriodAfter 1.53e-12 ***
## StationCodeSTTD:FlowActionPeriodAfter 5.89e-06 ***
## StationCodeLIB:FlowActionPeriodAfter < 2e-16 ***
## StationCodeRVB:FlowActionPeriodAfter 1.27e-14 ***
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.621005
## StationCodeI80:Flow:FlowActionPeriodDuring 0.139151
## StationCodeLIS:Flow:FlowActionPeriodDuring 0.000271 ***
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.523617
## StationCodeLIB:Flow:FlowActionPeriodDuring 0.073440 .
## StationCodeRVB:Flow:FlowActionPeriodDuring 0.036769 *
## StationCodeRD22:Flow:FlowActionPeriodAfter 0.453163
## StationCodeI80:Flow:FlowActionPeriodAfter 0.031634 *
## StationCodeLIS:Flow:FlowActionPeriodAfter < 2e-16 ***
## StationCodeSTTD:Flow:FlowActionPeriodAfter 3.90e-13 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter 1.64e-11 ***
## StationCodeRVB:Flow:FlowActionPeriodAfter 2.28e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.863 3.989 47.56 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.88 Deviance explained = 88.2%
## -REML = 1476.1 Scale est. = 0.11868 n = 3478appraise(m_gam_q_fa)
k.check(m_gam_q_fa)## k' edf k-index p-value
## s(DOY) 4 3.862781 0.9687532 0.0375draw(m_gam_q_fa, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_q_fa, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_q_fa))
Box.test(residuals(m_gam_q_fa), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_q_fa)
## X-squared = 7116.1, df = 20, p-value < 2.2e-16Model with AR() correlation structure
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_q_fa <- as.formula("Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + s(DOY, k = 5)")
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q_fa <- gamm(
f_gam_q_fa,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_q_fa_ar1 <- gamm(
f_gam_q_fa,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_fa_ar2 <- gamm(
f_gam_q_fa,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_fa_ar3 <- gamm(
f_gam_q_fa,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_q_fa$lme,
m_gamm_q_fa_ar1$lme,
m_gamm_q_fa_ar2$lme,
m_gamm_q_fa_ar3$lme
)## Model df AIC BIC logLik Test L.Ratio
## m_gamm_q_fa$lme 1 73 3098.134 3545.886 -1476.067
## m_gamm_q_fa_ar1$lme 2 74 -1985.717 -1531.832 1066.859 1 vs 2 5085.851
## m_gamm_q_fa_ar2$lme 3 75 -2005.949 -1545.930 1077.975 2 vs 3 22.232
## m_gamm_q_fa_ar3$lme 4 76 -2004.214 -1538.061 1078.107 3 vs 4 0.265
## p-value
## m_gamm_q_fa$lme
## m_gamm_q_fa_ar1$lme <.0001
## m_gamm_q_fa_ar2$lme <.0001
## m_gamm_q_fa_ar3$lme 0.6067Again, it looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
AR(2) Model
# Pull out the residuals and the GAM model
resid_q_fa_ar2 <- residuals(m_gamm_q_fa_ar2$lme, type = "normalized")
m_gamm_q_fa_ar2_gam <- m_gamm_q_fa_ar2$gam
summary(m_gamm_q_fa_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value
## (Intercept) 9.140e+00 3.300e-01 27.694
## Year_fct2016 -8.591e-02 4.210e-01 -0.204
## Year_fct2017 3.071e-01 3.994e-01 0.769
## Year_fct2018 1.029e-01 3.897e-01 0.264
## Year_fct2019 6.427e-02 3.912e-01 0.164
## StationCodeRD22 2.084e-01 4.107e-01 0.508
## StationCodeI80 3.448e-01 4.231e-01 0.815
## StationCodeLIS -6.571e-01 3.993e-01 -1.646
## StationCodeSTTD -2.267e-01 4.099e-01 -0.553
## StationCodeLIB -9.494e-01 4.148e-01 -2.289
## StationCodeRVB -1.386e+00 4.029e-01 -3.441
## Flow -1.686e-04 3.656e-04 -0.461
## FlowActionPeriodDuring 8.981e-02 1.096e-01 0.819
## FlowActionPeriodAfter 1.447e-01 1.313e-01 1.102
## Year_fct2016:StationCodeRD22 -4.345e-01 5.402e-01 -0.804
## Year_fct2017:StationCodeRD22 -5.513e-01 5.127e-01 -1.075
## Year_fct2018:StationCodeRD22 -3.939e-01 5.015e-01 -0.785
## Year_fct2019:StationCodeRD22 -1.286e-01 5.039e-01 -0.255
## Year_fct2016:StationCodeI80 1.236e-01 5.486e-01 0.225
## Year_fct2017:StationCodeI80 -3.138e-01 5.212e-01 -0.602
## Year_fct2018:StationCodeI80 -5.976e-01 5.102e-01 -1.171
## Year_fct2019:StationCodeI80 4.697e-01 5.125e-01 0.916
## Year_fct2016:StationCodeLIS 7.761e-01 5.260e-01 1.475
## Year_fct2017:StationCodeLIS 2.046e-01 5.104e-01 0.401
## Year_fct2018:StationCodeLIS -1.434e-01 4.971e-01 -0.288
## Year_fct2019:StationCodeLIS 4.044e-01 5.029e-01 0.804
## Year_fct2016:StationCodeSTTD 3.971e-01 5.426e-01 0.732
## Year_fct2017:StationCodeSTTD -6.679e-01 5.490e-01 -1.217
## Year_fct2018:StationCodeSTTD -9.501e-01 5.208e-01 -1.824
## Year_fct2019:StationCodeSTTD -1.144e+00 5.126e-01 -2.232
## Year_fct2016:StationCodeLIB 9.970e-01 5.295e-01 1.883
## Year_fct2017:StationCodeLIB -4.086e-01 5.154e-01 -0.793
## Year_fct2018:StationCodeLIB -2.176e+00 5.607e-01 -3.881
## Year_fct2019:StationCodeLIB -9.519e-01 5.694e-01 -1.672
## Year_fct2016:StationCodeRVB 3.678e-01 5.390e-01 0.682
## Year_fct2017:StationCodeRVB -5.412e-01 5.064e-01 -1.069
## Year_fct2018:StationCodeRVB -4.418e-01 4.948e-01 -0.893
## Year_fct2019:StationCodeRVB -8.452e-01 5.035e-01 -1.678
## StationCodeRD22:Flow 6.557e-04 1.189e-03 0.552
## StationCodeI80:Flow 4.714e-03 1.602e-03 2.943
## StationCodeLIS:Flow -8.155e-04 8.945e-04 -0.912
## StationCodeSTTD:Flow -1.186e-04 1.026e-03 -0.116
## StationCodeLIB:Flow 4.674e-04 3.725e-04 1.255
## StationCodeRVB:Flow 1.588e-04 3.657e-04 0.434
## Flow:FlowActionPeriodDuring -5.381e-05 3.758e-04 -0.143
## Flow:FlowActionPeriodAfter 7.104e-05 7.263e-04 0.098
## StationCodeRD22:FlowActionPeriodDuring 3.777e-02 1.509e-01 0.250
## StationCodeI80:FlowActionPeriodDuring 3.444e-01 1.540e-01 2.236
## StationCodeLIS:FlowActionPeriodDuring -1.107e-01 1.328e-01 -0.834
## StationCodeSTTD:FlowActionPeriodDuring 7.934e-03 1.340e-01 0.059
## StationCodeLIB:FlowActionPeriodDuring -4.421e-01 1.783e-01 -2.480
## StationCodeRVB:FlowActionPeriodDuring -1.782e-01 1.846e-01 -0.965
## StationCodeRD22:FlowActionPeriodAfter 8.258e-02 1.954e-01 0.423
## StationCodeI80:FlowActionPeriodAfter 4.533e-01 1.984e-01 2.285
## StationCodeLIS:FlowActionPeriodAfter 1.976e-02 1.635e-01 0.121
## StationCodeSTTD:FlowActionPeriodAfter -4.926e-01 1.658e-01 -2.972
## StationCodeLIB:FlowActionPeriodAfter -4.578e-01 1.974e-01 -2.320
## StationCodeRVB:FlowActionPeriodAfter -2.304e-01 1.999e-01 -1.153
## StationCodeRD22:Flow:FlowActionPeriodDuring -1.040e-03 1.172e-03 -0.887
## StationCodeI80:Flow:FlowActionPeriodDuring -5.914e-03 1.569e-03 -3.771
## StationCodeLIS:Flow:FlowActionPeriodDuring 1.118e-03 9.342e-04 1.197
## StationCodeSTTD:Flow:FlowActionPeriodDuring 9.933e-04 1.061e-03 0.936
## StationCodeLIB:Flow:FlowActionPeriodDuring -1.879e-04 3.889e-04 -0.483
## StationCodeRVB:Flow:FlowActionPeriodDuring 6.217e-05 3.761e-04 0.165
## StationCodeRD22:Flow:FlowActionPeriodAfter -4.517e-05 2.549e-03 -0.018
## StationCodeI80:Flow:FlowActionPeriodAfter -8.338e-03 2.762e-03 -3.019
## StationCodeLIS:Flow:FlowActionPeriodAfter -4.819e-05 1.269e-03 -0.038
## StationCodeSTTD:Flow:FlowActionPeriodAfter 4.786e-03 1.381e-03 3.465
## StationCodeLIB:Flow:FlowActionPeriodAfter -2.209e-04 7.319e-04 -0.302
## StationCodeRVB:Flow:FlowActionPeriodAfter -6.675e-05 7.263e-04 -0.092
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year_fct2016 0.838299
## Year_fct2017 0.441995
## Year_fct2018 0.791690
## Year_fct2019 0.869520
## StationCodeRD22 0.611816
## StationCodeI80 0.415252
## StationCodeLIS 0.099916 .
## StationCodeSTTD 0.580352
## StationCodeLIB 0.022136 *
## StationCodeRVB 0.000587 ***
## Flow 0.644739
## FlowActionPeriodDuring 0.412618
## FlowActionPeriodAfter 0.270588
## Year_fct2016:StationCodeRD22 0.421258
## Year_fct2017:StationCodeRD22 0.282374
## Year_fct2018:StationCodeRD22 0.432237
## Year_fct2019:StationCodeRD22 0.798634
## Year_fct2016:StationCodeI80 0.821794
## Year_fct2017:StationCodeI80 0.547132
## Year_fct2018:StationCodeI80 0.241558
## Year_fct2019:StationCodeI80 0.359519
## Year_fct2016:StationCodeLIS 0.140224
## Year_fct2017:StationCodeLIS 0.688545
## Year_fct2018:StationCodeLIS 0.772996
## Year_fct2019:StationCodeLIS 0.421423
## Year_fct2016:StationCodeSTTD 0.464298
## Year_fct2017:StationCodeSTTD 0.223845
## Year_fct2018:StationCodeSTTD 0.068178 .
## Year_fct2019:StationCodeSTTD 0.025692 *
## Year_fct2016:StationCodeLIB 0.059765 .
## Year_fct2017:StationCodeLIB 0.427980
## Year_fct2018:StationCodeLIB 0.000106 ***
## Year_fct2019:StationCodeLIB 0.094640 .
## Year_fct2016:StationCodeRVB 0.495041
## Year_fct2017:StationCodeRVB 0.285267
## Year_fct2018:StationCodeRVB 0.371985
## Year_fct2019:StationCodeRVB 0.093350 .
## StationCodeRD22:Flow 0.581215
## StationCodeI80:Flow 0.003268 **
## StationCodeLIS:Flow 0.362050
## StationCodeSTTD:Flow 0.908046
## StationCodeLIB:Flow 0.209715
## StationCodeRVB:Flow 0.664224
## Flow:FlowActionPeriodDuring 0.886142
## Flow:FlowActionPeriodAfter 0.922079
## StationCodeRD22:FlowActionPeriodDuring 0.802365
## StationCodeI80:FlowActionPeriodDuring 0.025409 *
## StationCodeLIS:FlowActionPeriodDuring 0.404406
## StationCodeSTTD:FlowActionPeriodDuring 0.952780
## StationCodeLIB:FlowActionPeriodDuring 0.013178 *
## StationCodeRVB:FlowActionPeriodDuring 0.334379
## StationCodeRD22:FlowActionPeriodAfter 0.672656
## StationCodeI80:FlowActionPeriodAfter 0.022380 *
## StationCodeLIS:FlowActionPeriodAfter 0.903816
## StationCodeSTTD:FlowActionPeriodAfter 0.002984 **
## StationCodeLIB:FlowActionPeriodAfter 0.020421 *
## StationCodeRVB:FlowActionPeriodAfter 0.249071
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.374929
## StationCodeI80:Flow:FlowActionPeriodDuring 0.000166 ***
## StationCodeLIS:Flow:FlowActionPeriodDuring 0.231434
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.349200
## StationCodeLIB:Flow:FlowActionPeriodDuring 0.629007
## StationCodeRVB:Flow:FlowActionPeriodDuring 0.868712
## StationCodeRD22:Flow:FlowActionPeriodAfter 0.985860
## StationCodeI80:Flow:FlowActionPeriodAfter 0.002555 **
## StationCodeLIS:Flow:FlowActionPeriodAfter 0.969721
## StationCodeSTTD:Flow:FlowActionPeriodAfter 0.000537 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter 0.762798
## StationCodeRVB:Flow:FlowActionPeriodAfter 0.926785
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.651 3.651 9.374 1.24e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.806
## Scale est. = 0.25102 n = 3478appraise(m_gamm_q_fa_ar2_gam)
k.check(m_gamm_q_fa_ar2_gam)## k' edf k-index p-value
## s(DOY) 4 3.650986 0.8149469 0draw(m_gamm_q_fa_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_q_fa_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_q_fa_ar2)
Box.test(resid_q_fa_ar2, lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: resid_q_fa_ar2
## X-squared = 30.35, df = 20, p-value = 0.06438The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_fa_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 0.342 0.849425
## StationCode 6 6.770 3.80e-07
## Flow 1 0.213 0.644739
## FlowActionPeriod 2 0.607 0.544928
## Year_fct:StationCode 24 3.414 3.87e-08
## StationCode:Flow 6 4.543 0.000135
## Flow:FlowActionPeriod 2 0.026 0.974524
## StationCode:FlowActionPeriod 12 5.610 1.29e-09
## StationCode:Flow:FlowActionPeriod 12 5.040 2.22e-08
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.651 3.651 9.374 1.24e-06Despite its complexity, we’ll use this model for the model selection process.
GAM Model 6: Flow with flow action period and water temperature
Initial Model
Finally, we’ll try running a GAM using the model structure from the flow with flow action period model and adding two-way interactions between Year:WaterTemp and Station:WaterTemp. First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_q_fa_wt <- gam(
Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)Lets look at the model summary and diagnostics:
summary(m_gam_q_fa_wt)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value
## (Intercept) 9.547e+00 2.926e-01 32.628
## Year_fct2016 -3.523e-01 3.056e-01 -1.153
## Year_fct2017 -5.764e-03 1.687e-01 -0.034
## Year_fct2018 -3.353e-01 1.608e-01 -2.086
## Year_fct2019 4.571e-01 1.558e-01 2.933
## StationCodeRD22 -4.980e-01 3.355e-01 -1.484
## StationCodeI80 7.357e-01 3.224e-01 2.281
## StationCodeLIS -2.400e+00 2.969e-01 -8.084
## StationCodeSTTD -1.276e+00 3.227e-01 -3.955
## StationCodeLIB -1.247e+00 4.001e-01 -3.117
## StationCodeRVB -2.148e+00 3.514e-01 -6.113
## Flow -7.271e-04 6.433e-04 -1.130
## FlowActionPeriodDuring 3.434e-02 9.046e-02 0.380
## FlowActionPeriodAfter 4.980e-01 1.057e-01 4.710
## WaterTemp -1.568e-02 1.146e-02 -1.368
## Year_fct2016:StationCodeRD22 -4.571e-01 1.005e-01 -4.548
## Year_fct2017:StationCodeRD22 -5.739e-01 8.460e-02 -6.783
## Year_fct2018:StationCodeRD22 -2.862e-01 8.045e-02 -3.558
## Year_fct2019:StationCodeRD22 -6.503e-02 8.367e-02 -0.777
## Year_fct2016:StationCodeI80 3.141e-01 1.040e-01 3.022
## Year_fct2017:StationCodeI80 -1.949e-01 8.623e-02 -2.260
## Year_fct2018:StationCodeI80 -3.712e-01 8.364e-02 -4.438
## Year_fct2019:StationCodeI80 6.428e-01 8.528e-02 7.538
## Year_fct2016:StationCodeLIS 1.408e-01 9.922e-02 1.419
## Year_fct2017:StationCodeLIS 1.859e-01 8.511e-02 2.185
## Year_fct2018:StationCodeLIS 1.916e-02 7.961e-02 0.241
## Year_fct2019:StationCodeLIS 6.092e-01 8.473e-02 7.190
## Year_fct2016:StationCodeSTTD 1.687e-01 1.025e-01 1.647
## Year_fct2017:StationCodeSTTD -5.588e-02 9.667e-02 -0.578
## Year_fct2018:StationCodeSTTD -5.559e-01 8.433e-02 -6.592
## Year_fct2019:StationCodeSTTD -1.018e+00 8.731e-02 -11.662
## Year_fct2016:StationCodeLIB 8.683e-01 1.027e-01 8.458
## Year_fct2017:StationCodeLIB -2.345e-01 9.240e-02 -2.538
## Year_fct2018:StationCodeLIB -1.994e+00 1.094e-01 -18.217
## Year_fct2019:StationCodeLIB -9.703e-01 1.172e-01 -8.279
## Year_fct2016:StationCodeRVB 1.142e-01 1.092e-01 1.046
## Year_fct2017:StationCodeRVB -3.724e-01 9.999e-02 -3.724
## Year_fct2018:StationCodeRVB -2.509e-01 8.781e-02 -2.857
## Year_fct2019:StationCodeRVB -6.336e-01 9.285e-02 -6.824
## StationCodeRD22:Flow 6.978e-04 1.383e-03 0.505
## StationCodeI80:Flow 4.347e-04 1.870e-03 0.232
## StationCodeLIS:Flow 5.075e-03 9.885e-04 5.134
## StationCodeSTTD:Flow 3.462e-03 1.099e-03 3.149
## StationCodeLIB:Flow 1.214e-03 6.495e-04 1.870
## StationCodeRVB:Flow 7.217e-04 6.433e-04 1.122
## Flow:FlowActionPeriodDuring -1.792e-04 6.558e-04 -0.273
## Flow:FlowActionPeriodAfter -6.344e-03 1.178e-03 -5.385
## StationCodeRD22:FlowActionPeriodDuring -1.234e-01 1.156e-01 -1.068
## StationCodeI80:FlowActionPeriodDuring -8.921e-02 1.197e-01 -0.746
## StationCodeLIS:FlowActionPeriodDuring 5.136e-01 1.132e-01 4.539
## StationCodeSTTD:FlowActionPeriodDuring 5.521e-01 1.237e-01 4.463
## StationCodeLIB:FlowActionPeriodDuring -1.717e-01 1.636e-01 -1.050
## StationCodeRVB:FlowActionPeriodDuring -1.081e-01 1.414e-01 -0.764
## StationCodeRD22:FlowActionPeriodAfter -2.807e-01 1.429e-01 -1.965
## StationCodeI80:FlowActionPeriodAfter -5.002e-01 1.462e-01 -3.422
## StationCodeLIS:FlowActionPeriodAfter -4.052e-02 1.220e-01 -0.332
## StationCodeSTTD:FlowActionPeriodAfter -1.510e-01 1.302e-01 -1.160
## StationCodeLIB:FlowActionPeriodAfter -1.600e+00 1.802e-01 -8.878
## StationCodeRVB:FlowActionPeriodAfter -6.906e-01 1.593e-01 -4.336
## Year_fct2016:WaterTemp 1.554e-02 1.201e-02 1.294
## Year_fct2017:WaterTemp 7.768e-03 6.761e-03 1.149
## Year_fct2018:WaterTemp 1.534e-02 6.559e-03 2.339
## Year_fct2019:WaterTemp -2.449e-02 6.218e-03 -3.938
## StationCodeRD22:WaterTemp 4.178e-02 1.249e-02 3.345
## StationCodeI80:WaterTemp -2.552e-03 1.187e-02 -0.215
## StationCodeLIS:WaterTemp 8.755e-02 1.119e-02 7.823
## StationCodeSTTD:WaterTemp 2.874e-02 1.227e-02 2.343
## StationCodeLIB:WaterTemp 2.504e-02 1.611e-02 1.555
## StationCodeRVB:WaterTemp 3.059e-02 1.375e-02 2.224
## StationCodeRD22:Flow:FlowActionPeriodDuring -1.012e-03 1.393e-03 -0.727
## StationCodeI80:Flow:FlowActionPeriodDuring -1.383e-03 1.874e-03 -0.738
## StationCodeLIS:Flow:FlowActionPeriodDuring -5.526e-03 1.007e-03 -5.491
## StationCodeSTTD:Flow:FlowActionPeriodDuring -1.822e-03 1.117e-03 -1.631
## StationCodeLIB:Flow:FlowActionPeriodDuring -1.014e-05 6.641e-04 -0.015
## StationCodeRVB:Flow:FlowActionPeriodDuring 1.844e-04 6.559e-04 0.281
## StationCodeRD22:Flow:FlowActionPeriodAfter -4.200e-03 2.534e-03 -1.657
## StationCodeI80:Flow:FlowActionPeriodAfter -1.946e-03 2.783e-03 -0.699
## StationCodeLIS:Flow:FlowActionPeriodAfter 8.734e-03 1.560e-03 5.598
## StationCodeSTTD:Flow:FlowActionPeriodAfter 1.013e-02 1.662e-03 6.094
## StationCodeLIB:Flow:FlowActionPeriodAfter 5.708e-03 1.184e-03 4.820
## StationCodeRVB:Flow:FlowActionPeriodAfter 6.351e-03 1.178e-03 5.389
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year_fct2016 0.249091
## Year_fct2017 0.972744
## Year_fct2018 0.037076 *
## Year_fct2019 0.003375 **
## StationCodeRD22 0.137856
## StationCodeI80 0.022584 *
## StationCodeLIS 8.60e-16 ***
## StationCodeSTTD 7.82e-05 ***
## StationCodeLIB 0.001842 **
## StationCodeRVB 1.09e-09 ***
## Flow 0.258469
## FlowActionPeriodDuring 0.704248
## FlowActionPeriodAfter 2.58e-06 ***
## WaterTemp 0.171506
## Year_fct2016:StationCodeRD22 5.60e-06 ***
## Year_fct2017:StationCodeRD22 1.38e-11 ***
## Year_fct2018:StationCodeRD22 0.000379 ***
## Year_fct2019:StationCodeRD22 0.437081
## Year_fct2016:StationCodeI80 0.002533 **
## Year_fct2017:StationCodeI80 0.023879 *
## Year_fct2018:StationCodeI80 9.35e-06 ***
## Year_fct2019:StationCodeI80 6.10e-14 ***
## Year_fct2016:StationCodeLIS 0.155939
## Year_fct2017:StationCodeLIS 0.028987 *
## Year_fct2018:StationCodeLIS 0.809791
## Year_fct2019:StationCodeLIS 7.93e-13 ***
## Year_fct2016:StationCodeSTTD 0.099692 .
## Year_fct2017:StationCodeSTTD 0.563289
## Year_fct2018:StationCodeSTTD 5.03e-11 ***
## Year_fct2019:StationCodeSTTD < 2e-16 ***
## Year_fct2016:StationCodeLIB < 2e-16 ***
## Year_fct2017:StationCodeLIB 0.011181 *
## Year_fct2018:StationCodeLIB < 2e-16 ***
## Year_fct2019:StationCodeLIB < 2e-16 ***
## Year_fct2016:StationCodeRVB 0.295425
## Year_fct2017:StationCodeRVB 0.000199 ***
## Year_fct2018:StationCodeRVB 0.004298 **
## Year_fct2019:StationCodeRVB 1.05e-11 ***
## StationCodeRD22:Flow 0.613905
## StationCodeI80:Flow 0.816208
## StationCodeLIS:Flow 2.99e-07 ***
## StationCodeSTTD:Flow 0.001652 **
## StationCodeLIB:Flow 0.061637 .
## StationCodeRVB:Flow 0.261987
## Flow:FlowActionPeriodDuring 0.784664
## Flow:FlowActionPeriodAfter 7.73e-08 ***
## StationCodeRD22:FlowActionPeriodDuring 0.285703
## StationCodeI80:FlowActionPeriodDuring 0.455987
## StationCodeLIS:FlowActionPeriodDuring 5.86e-06 ***
## StationCodeSTTD:FlowActionPeriodDuring 8.33e-06 ***
## StationCodeLIB:FlowActionPeriodDuring 0.293788
## StationCodeRVB:FlowActionPeriodDuring 0.444642
## StationCodeRD22:FlowActionPeriodAfter 0.049498 *
## StationCodeI80:FlowActionPeriodAfter 0.000629 ***
## StationCodeLIS:FlowActionPeriodAfter 0.739835
## StationCodeSTTD:FlowActionPeriodAfter 0.246237
## StationCodeLIB:FlowActionPeriodAfter < 2e-16 ***
## StationCodeRVB:FlowActionPeriodAfter 1.49e-05 ***
## Year_fct2016:WaterTemp 0.195712
## Year_fct2017:WaterTemp 0.250655
## Year_fct2018:WaterTemp 0.019415 *
## Year_fct2019:WaterTemp 8.39e-05 ***
## StationCodeRD22:WaterTemp 0.000832 ***
## StationCodeI80:WaterTemp 0.829731
## StationCodeLIS:WaterTemp 6.85e-15 ***
## StationCodeSTTD:WaterTemp 0.019208 *
## StationCodeLIB:WaterTemp 0.120121
## StationCodeRVB:WaterTemp 0.026229 *
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.467471
## StationCodeI80:Flow:FlowActionPeriodDuring 0.460348
## StationCodeLIS:Flow:FlowActionPeriodDuring 4.30e-08 ***
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.103086
## StationCodeLIB:Flow:FlowActionPeriodDuring 0.987816
## StationCodeRVB:Flow:FlowActionPeriodDuring 0.778648
## StationCodeRD22:Flow:FlowActionPeriodAfter 0.097598 .
## StationCodeI80:Flow:FlowActionPeriodAfter 0.484364
## StationCodeLIS:Flow:FlowActionPeriodAfter 2.34e-08 ***
## StationCodeSTTD:Flow:FlowActionPeriodAfter 1.22e-09 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter 1.50e-06 ***
## StationCodeRVB:Flow:FlowActionPeriodAfter 7.55e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.848 3.988 38.09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.885 Deviance explained = 88.8%
## -REML = 1434.9 Scale est. = 0.11336 n = 3478appraise(m_gam_q_fa_wt)
k.check(m_gam_q_fa_wt)## k' edf k-index p-value
## s(DOY) 4 3.84787 0.9605145 0.015draw(m_gam_q_fa_wt, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_q_fa_wt, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_q_fa_wt))
Box.test(residuals(m_gam_q_fa_wt), lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: residuals(m_gam_q_fa_wt)
## X-squared = 6772.7, df = 20, p-value < 2.2e-16Model with AR() correlation structure
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_q_fa_wt <- as.formula(
"Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)"
)
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q_fa_wt <- gamm(
f_gam_q_fa_wt,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_q_fa_wt_ar1 <- gamm(
f_gam_q_fa_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_fa_wt_ar2 <- gamm(
f_gam_q_fa_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_fa_wt_ar3 <- gamm(
f_gam_q_fa_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_q_fa_wt$lme,
m_gamm_q_fa_wt_ar1$lme,
m_gamm_q_fa_wt_ar2$lme,
m_gamm_q_fa_wt_ar3$lme
)## Model df AIC BIC logLik Test L.Ratio
## m_gamm_q_fa_wt$lme 1 84 3037.848 3552.798 -1434.924
## m_gamm_q_fa_wt_ar1$lme 2 85 -1962.244 -1441.163 1066.122 1 vs 2 5002.092
## m_gamm_q_fa_wt_ar2$lme 3 86 -1985.590 -1458.380 1078.795 2 vs 3 25.347
## m_gamm_q_fa_wt_ar3$lme 4 87 -1983.981 -1450.641 1078.991 3 vs 4 0.391
## p-value
## m_gamm_q_fa_wt$lme
## m_gamm_q_fa_wt_ar1$lme <.0001
## m_gamm_q_fa_wt_ar2$lme <.0001
## m_gamm_q_fa_wt_ar3$lme 0.5318Again, it looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
AR(2) Model
# Pull out the residuals and the GAM model
resid_q_fa_wt_ar2 <- residuals(m_gamm_q_fa_wt_ar2$lme, type = "normalized")
m_gamm_q_fa_wt_ar2_gam <- m_gamm_q_fa_wt_ar2$gam
summary(m_gamm_q_fa_wt_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value
## (Intercept) 8.796e+00 6.363e-01 13.824
## Year_fct2016 1.238e+00 6.561e-01 1.887
## Year_fct2017 1.791e-01 5.995e-01 0.299
## Year_fct2018 -4.588e-01 6.028e-01 -0.761
## Year_fct2019 -1.297e-01 5.750e-01 -0.226
## StationCodeRD22 -1.322e+00 6.714e-01 -1.969
## StationCodeI80 4.878e-01 6.551e-01 0.745
## StationCodeLIS -1.669e+00 6.586e-01 -2.534
## StationCodeSTTD -3.777e-01 6.810e-01 -0.555
## StationCodeLIB -3.107e-01 7.561e-01 -0.411
## StationCodeRVB -1.787e+00 8.259e-01 -2.164
## Flow -2.257e-04 3.659e-04 -0.617
## FlowActionPeriodDuring 7.161e-02 1.103e-01 0.649
## FlowActionPeriodAfter 1.159e-01 1.356e-01 0.855
## WaterTemp 1.664e-02 2.240e-02 0.743
## Year_fct2016:StationCodeRD22 -6.248e-01 5.824e-01 -1.073
## Year_fct2017:StationCodeRD22 -6.087e-01 5.523e-01 -1.102
## Year_fct2018:StationCodeRD22 -3.091e-01 5.418e-01 -0.571
## Year_fct2019:StationCodeRD22 -1.262e-01 5.437e-01 -0.232
## Year_fct2016:StationCodeI80 1.584e-01 5.929e-01 0.267
## Year_fct2017:StationCodeI80 -2.583e-01 5.628e-01 -0.459
## Year_fct2018:StationCodeI80 -5.353e-01 5.522e-01 -0.969
## Year_fct2019:StationCodeI80 4.972e-01 5.542e-01 0.897
## Year_fct2016:StationCodeLIS 6.262e-01 5.677e-01 1.103
## Year_fct2017:StationCodeLIS 1.909e-01 5.499e-01 0.347
## Year_fct2018:StationCodeLIS -1.089e-01 5.368e-01 -0.203
## Year_fct2019:StationCodeLIS 3.895e-01 5.428e-01 0.718
## Year_fct2016:StationCodeSTTD 3.693e-01 5.864e-01 0.630
## Year_fct2017:StationCodeSTTD -6.079e-01 5.920e-01 -1.027
## Year_fct2018:StationCodeSTTD -9.043e-01 5.621e-01 -1.609
## Year_fct2019:StationCodeSTTD -1.126e+00 5.536e-01 -2.034
## Year_fct2016:StationCodeLIB 9.232e-01 5.719e-01 1.614
## Year_fct2017:StationCodeLIB -3.232e-01 5.559e-01 -0.581
## Year_fct2018:StationCodeLIB -2.101e+00 6.040e-01 -3.478
## Year_fct2019:StationCodeLIB -8.633e-01 6.117e-01 -1.411
## Year_fct2016:StationCodeRVB 2.471e-01 5.819e-01 0.425
## Year_fct2017:StationCodeRVB -4.778e-01 5.461e-01 -0.875
## Year_fct2018:StationCodeRVB -3.624e-01 5.353e-01 -0.677
## Year_fct2019:StationCodeRVB -7.749e-01 5.446e-01 -1.423
## StationCodeRD22:Flow 6.934e-04 1.181e-03 0.587
## StationCodeI80:Flow 4.800e-03 1.591e-03 3.017
## StationCodeLIS:Flow -2.042e-04 9.037e-04 -0.226
## StationCodeSTTD:Flow 1.321e-04 1.033e-03 0.128
## StationCodeLIB:Flow 4.830e-04 3.730e-04 1.295
## StationCodeRVB:Flow 2.179e-04 3.660e-04 0.595
## Flow:FlowActionPeriodDuring 2.747e-05 3.769e-04 0.073
## Flow:FlowActionPeriodAfter 6.981e-05 7.286e-04 0.096
## StationCodeRD22:FlowActionPeriodDuring 7.152e-02 1.515e-01 0.472
## StationCodeI80:FlowActionPeriodDuring 3.482e-01 1.549e-01 2.248
## StationCodeLIS:FlowActionPeriodDuring -7.908e-02 1.334e-01 -0.593
## StationCodeSTTD:FlowActionPeriodDuring 9.413e-03 1.343e-01 0.070
## StationCodeLIB:FlowActionPeriodDuring -4.343e-01 1.791e-01 -2.426
## StationCodeRVB:FlowActionPeriodDuring -1.575e-01 1.846e-01 -0.854
## StationCodeRD22:FlowActionPeriodAfter 2.049e-01 2.013e-01 1.018
## StationCodeI80:FlowActionPeriodAfter 4.531e-01 2.040e-01 2.221
## StationCodeLIS:FlowActionPeriodAfter 1.002e-01 1.688e-01 0.594
## StationCodeSTTD:FlowActionPeriodAfter -4.901e-01 1.705e-01 -2.874
## StationCodeLIB:FlowActionPeriodAfter -4.527e-01 2.014e-01 -2.247
## StationCodeRVB:FlowActionPeriodAfter -1.957e-01 2.062e-01 -0.949
## Year_fct2016:WaterTemp -5.317e-02 2.012e-02 -2.642
## Year_fct2017:WaterTemp 4.584e-03 1.825e-02 0.251
## Year_fct2018:WaterTemp 2.569e-02 1.940e-02 1.325
## Year_fct2019:WaterTemp 8.629e-03 1.730e-02 0.499
## StationCodeRD22:WaterTemp 6.962e-02 2.193e-02 3.175
## StationCodeI80:WaterTemp -8.944e-03 2.043e-02 -0.438
## StationCodeLIS:WaterTemp 4.623e-02 2.165e-02 2.135
## StationCodeSTTD:WaterTemp 5.668e-03 2.295e-02 0.247
## StationCodeLIB:WaterTemp -3.515e-02 2.818e-02 -1.248
## StationCodeRVB:WaterTemp 1.683e-02 3.139e-02 0.536
## StationCodeRD22:Flow:FlowActionPeriodDuring -1.121e-03 1.164e-03 -0.963
## StationCodeI80:Flow:FlowActionPeriodDuring -5.969e-03 1.557e-03 -3.832
## StationCodeLIS:Flow:FlowActionPeriodDuring 4.993e-04 9.440e-04 0.529
## StationCodeSTTD:Flow:FlowActionPeriodDuring 7.451e-04 1.066e-03 0.699
## StationCodeLIB:Flow:FlowActionPeriodDuring -2.590e-04 3.895e-04 -0.665
## StationCodeRVB:Flow:FlowActionPeriodDuring -1.933e-05 3.772e-04 -0.051
## StationCodeRD22:Flow:FlowActionPeriodAfter -9.655e-04 2.553e-03 -0.378
## StationCodeI80:Flow:FlowActionPeriodAfter -7.980e-03 2.766e-03 -2.885
## StationCodeLIS:Flow:FlowActionPeriodAfter -6.062e-04 1.275e-03 -0.475
## StationCodeSTTD:Flow:FlowActionPeriodAfter 4.649e-03 1.391e-03 3.343
## StationCodeLIB:Flow:FlowActionPeriodAfter -1.939e-04 7.341e-04 -0.264
## StationCodeRVB:Flow:FlowActionPeriodAfter -6.629e-05 7.287e-04 -0.091
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year_fct2016 0.059184 .
## Year_fct2017 0.765131
## Year_fct2018 0.446654
## Year_fct2019 0.821571
## StationCodeRD22 0.049085 *
## StationCodeI80 0.456585
## StationCodeLIS 0.011330 *
## StationCodeSTTD 0.579142
## StationCodeLIB 0.681142
## StationCodeRVB 0.030529 *
## Flow 0.537306
## FlowActionPeriodDuring 0.516207
## FlowActionPeriodAfter 0.392757
## WaterTemp 0.457427
## Year_fct2016:StationCodeRD22 0.283391
## Year_fct2017:StationCodeRD22 0.270529
## Year_fct2018:StationCodeRD22 0.568362
## Year_fct2019:StationCodeRD22 0.816454
## Year_fct2016:StationCodeI80 0.789413
## Year_fct2017:StationCodeI80 0.646245
## Year_fct2018:StationCodeI80 0.332393
## Year_fct2019:StationCodeI80 0.369731
## Year_fct2016:StationCodeLIS 0.270033
## Year_fct2017:StationCodeLIS 0.728490
## Year_fct2018:StationCodeLIS 0.839244
## Year_fct2019:StationCodeLIS 0.473020
## Year_fct2016:StationCodeSTTD 0.528925
## Year_fct2017:StationCodeSTTD 0.304603
## Year_fct2018:StationCodeSTTD 0.107788
## Year_fct2019:StationCodeSTTD 0.042061 *
## Year_fct2016:StationCodeLIB 0.106552
## Year_fct2017:StationCodeLIB 0.560996
## Year_fct2018:StationCodeLIB 0.000511 ***
## Year_fct2019:StationCodeLIB 0.158224
## Year_fct2016:StationCodeRVB 0.671099
## Year_fct2017:StationCodeRVB 0.381721
## Year_fct2018:StationCodeRVB 0.498377
## Year_fct2019:StationCodeRVB 0.154848
## StationCodeRD22:Flow 0.557060
## StationCodeI80:Flow 0.002576 **
## StationCodeLIS:Flow 0.821271
## StationCodeSTTD:Flow 0.898262
## StationCodeLIB:Flow 0.195464
## StationCodeRVB:Flow 0.551666
## Flow:FlowActionPeriodDuring 0.941913
## Flow:FlowActionPeriodAfter 0.923674
## StationCodeRD22:FlowActionPeriodDuring 0.636993
## StationCodeI80:FlowActionPeriodDuring 0.024636 *
## StationCodeLIS:FlowActionPeriodDuring 0.553410
## StationCodeSTTD:FlowActionPeriodDuring 0.944114
## StationCodeLIB:FlowActionPeriodDuring 0.015339 *
## StationCodeRVB:FlowActionPeriodDuring 0.393441
## StationCodeRD22:FlowActionPeriodAfter 0.308635
## StationCodeI80:FlowActionPeriodAfter 0.026442 *
## StationCodeLIS:FlowActionPeriodAfter 0.552818
## StationCodeSTTD:FlowActionPeriodAfter 0.004072 **
## StationCodeLIB:FlowActionPeriodAfter 0.024681 *
## StationCodeRVB:FlowActionPeriodAfter 0.342716
## Year_fct2016:WaterTemp 0.008277 **
## Year_fct2017:WaterTemp 0.801701
## Year_fct2018:WaterTemp 0.185344
## Year_fct2019:WaterTemp 0.618037
## StationCodeRD22:WaterTemp 0.001511 **
## StationCodeI80:WaterTemp 0.661513
## StationCodeLIS:WaterTemp 0.032853 *
## StationCodeSTTD:WaterTemp 0.804969
## StationCodeLIB:WaterTemp 0.212274
## StationCodeRVB:WaterTemp 0.591847
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.335423
## StationCodeI80:Flow:FlowActionPeriodDuring 0.000129 ***
## StationCodeLIS:Flow:FlowActionPeriodDuring 0.596880
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.484804
## StationCodeLIB:Flow:FlowActionPeriodDuring 0.506064
## StationCodeRVB:Flow:FlowActionPeriodDuring 0.959130
## StationCodeRD22:Flow:FlowActionPeriodAfter 0.705263
## StationCodeI80:Flow:FlowActionPeriodAfter 0.003938 **
## StationCodeLIS:Flow:FlowActionPeriodAfter 0.634486
## StationCodeSTTD:Flow:FlowActionPeriodAfter 0.000838 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter 0.791721
## StationCodeRVB:Flow:FlowActionPeriodAfter 0.927530
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.623 3.623 7.084 2.26e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.797
## Scale est. = 0.27259 n = 3478appraise(m_gamm_q_fa_wt_ar2_gam)
k.check(m_gamm_q_fa_wt_ar2_gam)## k' edf k-index p-value
## s(DOY) 4 3.622926 0.8021592 0draw(m_gamm_q_fa_wt_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_q_fa_wt_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_q_fa_wt_ar2)
Box.test(resid_q_fa_wt_ar2, lag = 20, type = 'Ljung-Box')##
## Box-Ljung test
##
## data: resid_q_fa_wt_ar2
## X-squared = 30.263, df = 20, p-value = 0.0657The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_fa_wt_ar2_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 2.098 0.078512
## StationCode 6 4.328 0.000234
## Flow 1 0.381 0.537306
## FlowActionPeriod 2 0.366 0.693855
## WaterTemp 1 0.552 0.457427
## Year_fct:StationCode 24 2.863 4.02e-06
## StationCode:Flow 6 3.598 0.001470
## Flow:FlowActionPeriod 2 0.005 0.994855
## StationCode:FlowActionPeriod 12 5.902 2.95e-10
## Year_fct:WaterTemp 4 4.121 0.002471
## StationCode:WaterTemp 6 5.271 2.04e-05
## StationCode:Flow:FlowActionPeriod 12 4.904 4.35e-08
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.623 3.623 7.084 2.26e-05Despite its complexity, we’ll use this model for the model selection process.
Model selection - GAM models
As a summary, here are the 6 GAM models we are comparing:
- Model 1 -
m_gamm_cat_ar2- Categorical predictors only
Formula: Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, k = 5) - Model 2 -
m_gamm_cat_wt_ar2b- Categorical predictors adding water temperature
Formula: Chla_log ~ (Year_fct + StationCode + WaterTemp)^2 + StationCode * FlowActionPeriod + s(DOY, k = 5) - Model 3 -
m_gamm_q_ar2b- Flow as continuous predictor
Formula: Chla_log ~ Year_fct * StationCode + Flow * StationCode + s(DOY, k = 5) - Model 4 -
m_gamm_q_wt_ar2- Flow and water temperature as continuous predictors
Formula: Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5) - Model 5 -
m_gamm_q_fa_ar2- Flow with flow action period
Formula: Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + s(DOY, k = 5) - Model 6 -
m_gamm_q_fa_wt_ar2- Flow with flow action period and water temperature
Formula: Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
# AIC values
df_gamm_aic <-
AIC(
m_gamm_cat_ar2$lme,
m_gamm_cat_wt_ar2b$lme,
m_gamm_q_ar2b$lme,
m_gamm_q_wt_ar2$lme,
m_gamm_q_fa_ar2$lme,
m_gamm_q_fa_wt_ar2$lme
) %>%
as_tibble(rownames = "gamm_model")
# BIC values
df_gamm_bic <-
BIC(
m_gamm_cat_ar2$lme,
m_gamm_cat_wt_ar2b$lme,
m_gamm_q_ar2b$lme,
m_gamm_q_wt_ar2$lme,
m_gamm_q_fa_ar2$lme,
m_gamm_q_fa_wt_ar2$lme
) %>%
as_tibble(rownames = "gamm_model")
# Combine AIC and BIC and display results
left_join(df_gamm_aic, df_gamm_bic, by = join_by(gamm_model, df))## # A tibble: 6 × 4
## gamm_model df AIC BIC
## <chr> <dbl> <dbl> <dbl>
## 1 m_gamm_cat_ar2$lme 62 -2189. -1809.
## 2 m_gamm_cat_wt_ar2b$lme 65 -2204. -1806.
## 3 m_gamm_q_ar2b$lme 47 -2198. -1909.
## 4 m_gamm_q_wt_ar2$lme 58 -2174. -1818.
## 5 m_gamm_q_fa_ar2$lme 75 -2006. -1546.
## 6 m_gamm_q_fa_wt_ar2$lme 86 -1986. -1458.According to AIC, model 2 (categorical predictors adding water temperature) was the best model, while BIC preferred model 3 (flow as continuous predictor). We’ll look at these two models more closely.
m_gamm_cat_wt_ar2b %>% saveRDS(here("manuscript_synthesis/model/interim/chla_cat_wt_gamm_model.rds"))
m_gamm_q_ar2b %>% saveRDS(here("manuscript_synthesis/model/interim/chla_flow_gamm_model.rds"))
df_chla_wt_q %>% saveRDS(here("manuscript_synthesis/model/interim/chla_wt_flow_model_data.rds"))